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The evolution of star clusters is studied using A'^body simulations in which the 
evolution of single stars and binaries are taken self-consistently into account. Initial 
conditions are chosen to represent relatively young Galactic open clusters, such as the 
■ Pleiades, Praesepe and the Hyades. The calculations include a realistic mass function, 

pinj primordial binaries and the external potential of the parent Galaxy. 

Our model clusters are generally significantly flattened in the Galactic tidal field, 
\Z^ ' and dissolve before deep core collapse occurs. The binary fraction decreases initially 

^ \ due to the destruction of soft binaries, but increases later because lower mass single 

. stars escape more easily than the more massive binaries. At late times, the cluster 

core is quite rich in giants and white dwarfs. There is no evidence for preferential 
evaporation of old white dwarfs, on the contrary the formed white dwarfs are likely 
to remain in the cluster. Stars tend to escape from the cluster through the first and 
5— i ' second Lagrange points, in the direction of and away from the Galactic center. 

Mass segregation manifests itself in our models well within an initial relaxation 
time. As expected, giants and white dwarfs are much more strongly affected by mass 
segregation than main-sequence stars. However, the clusters to which we compare our 
results all appear to be somewhat more affected by mass segregation than our models 
would suggest. 

Open clusters are dynamically rather inactive. However, the combined effect of 
stellar mass loss and evaporation of stars from the cluster potential drives its disso- 
lution on a much shorter timescale than if these effects are neglected. The often-used 
argument that a star cluster is barely older than its relaxation time and therefore 
cannot be dynamically evolved is clearly in error for the majority of star clusters. 

An observation of a blue straggler in an eccentric orbit around an unevolved star 
or a blue straggler of more than twice the turn off mass might indicate past dynamical 
activity. We find two distinct populations of blue stragglers: those formed above the 
main-sequence turn off, and those which appear as blue stragglers as the cluster's 
turnoff drops below the mass of the rejuvenated star. 

1 INTRODUCTION 

Star clusters are remarkable laboratories for the study of 
fundamental astrophysical processes spanning the range of 
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stellar evolution, from birth to death. Open clusters are 
relatively small (typically ^ 10* stars), young (generally 
^ 1 Gyr old) systems found primarily in the Galactic disk. 
Roughly 1100 are known to lie within a few kiloparsecs of the 
Sun (Lynga 1987a). As a class, they afford us the opportu- 
nity to witness star formation and explore the complex evo- 
lution of young stellar systems. Since stars are born in stax 
clusters and become part of the general Galactic population 
only when the parent cluster dissolves in the Galactic tidal 
field, it is of considerable interest to understand how this 
process occurs, on what time scale it operates, and what are 
the detailed observable consequences of cluster evolution. 

At least as a first approximation, many star clusters are 
quite well characterized dynamically as pure iV-body sys- 
tems, in that there are few collisions or close encounters 
between stars and the clusters themselves are relatively iso- 
lated in space. The same argument can be made for the evo- 
lution of the stars and binaries without interactions, which 
may also quite well explain the properties of open clusters. 
This simple picture rapidly becomes inadequate when both 
effects are combined — cluster evolution is actuality an intri- 
cate mix of dynamics, stellar evolution, and external (tidal) 
influences, and the subtle interplay between stellar dynamics 
and stellar physics makes for a formidable modeling prob- 
lem. Of the many physical processes influencing cluster evo- 
lution, probably the most important axe the effects of stellar 
evolution, the tidal field of the Galaxy, and the evolution and 
dynamics of binary stars. We present here a scries of mod- 
els in which all these effects are taken self-consistently into 
account. 

The processes just mentioned are tightly coupled, com- 
plicating the evolution and making it difficult to isolate the 

importance of each individual effect. We refer to Mcylan & 
Heggie (1997) for a recent review. Generally, speaking, we 
can say mass loss from stellar evolution is of greatest impor- 
tance during the first few tens of millions of years of cluster 
evolution, and may well result in the disruption of the entire 
cluster (cf. Chernoff & Weinberg 1990, Fukushige & Heggie 
1995, Takahashi & Portegies Zwart 2000). If the cluster sur- 
vives this early phase, stellar evolutionary time scales soon 
become longer than the time scales for dynamical evolution, 
and two-body relaxation and tidal effects become dominant. 
Ultimately, these effects cause the cluster to dissociate and 
the stars to become part of the general "field" population of 
the disk. 

There is strong observational evidence that star clus- 
ters contain substantial binary populations, quite possibly 
as rich as those found in the field (see Rubenstein 1997). 
The properties and numbers of observed cluster binaries 
cannot be explained by internal formation processes, such 
as three body dynamics or two body tidal capture (Aarseth 
& Lecar 1975). The majority must be primordial, i.e. formed 
together with the single stars at the cluster's birth. These 
observations axe important to cluster dynamics, as a clus- 
ter's evolution depends strongly on its binary population 
and even a small initial binary fraction can play a pivotal 
role in governing cluster dynamics (Goodman & Hut 1989, 



McMillan et al. 1990, 1991a and 1991b Gao et al. 1991, Hut 
et al. 1992, McMillan & Hut 1994). Binaries are also crucial 
to cluster stellar evolution. The possibility of mass transfer 
between binary components permits wholly new stellar evo- 
lutionary states to arise; in addition, the presence of binaries 
will enhance the rate of stellar collisions and close encoun- 
ters, through the temporary capture of single stars and other 
binaries in three-body resonance encounters (Verbunt & Hut 
1987; Portegies Zwart et al. 1998). 

Until quite recently, dynamical models have tended to 
exclude binaries, for the good practical reasons that (1) bi- 
naries slow down the calculations dramatically and tend to 
induce numerical errors, and (2) their internal evolution is 
much more complicated than the evolution of single stars. 
However, it has also long boon known that cluster models 
lacking adequate treatments of binary systems are of at best 
limited validity. In this paper we begin to address these lim- 
itations by performing calculations of open clusters contain- 
ing substantial numbers of primordial binaries, with the goal 
of studying the mutual influence of stellar evolution and stel- 
lar dynamics in these systems — the "ecology" of star clusters 
(Heggie 1992). 

The first paper in this series (Portegies Zwart et 
al. 1997a, hereafter paper I) quantified the effect of collisions 
on stellar evolution and attempted to assess the correspond- 
ing changes in the stellar population. The stellar number 
density was held constant in these calculations, thus exclud- 
ing the possibility of any interplay between the dynamical 
evolution of the cluster and collisions between stars. In the 
second paper in this series (Portegies Zwart et al. 1997b, 
hereafter paper II) , the evolution of a population of primor- 
dial binaries was followed in time by tracking in detail the 
results of encounters between single stars and binaries. The 
assumption of constant stellar number density was relaxed 
in Portegies Zwart et al. (1999, hereafter paper HI), where 
the dynamical evolution of a star cluster (without primordial 
binaries) was followed in detail using A^-body calculations. 

This paper continues the process of relaxing the simpli- 
fying assumptions made in Paper I. As in Paper HI, we cal- 
culate the dynamical evolution of our model system by direct 
(A''-body) integration of the system, but now including both 
stars and binaries in the computational mix. Binary evolu- 
tion is incorporated into the A/'-body treatment, accounting 
for changes in binary orbital parameters due to stellar mass 
loss, supernovae, tidal forces between the stars, mass trans- 
fer from one star to its companion, general relativistic cor- 
rections, etc. Encounters between binaries and single stars 
and higher order encounters (between binaries and binaries 
and between binaries and triples etc.) axe fully integrated 
as are the orbits of the other stars in the AT-body system. 
All changes in the stellar and binary population caused by 
stellar evolution are fed back into the dynamical evolution of 
the parent cluster, allowing stellar dynamics and the stellar 
evolution to be studied simultaneously and self-consistently. 

As an initial case study, we concentrate here on open 
star clusters near the Sun (between ~ 6 and lOkpc from the 
Galactic center), which are less than 1 billion years old, and 
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which initially contain a few thousand (~ 2000-3000) stars, 
roughly half of them in binaries. Well known clusters fit- 
ting this general description are the Pleiades, Praesepe and 
the Hyades. Such systems are small enough that multiple 
simulations can bo performed in order to improve statistical 
coverage of their properties, yet they are large enough and 
old enough that both stellar evolution and stellar dynamics 
have had time to play significant roles in determining their 
present structure and appearance. 

In an effort to close the gap between theoretical and 
observational studies of cluster structure, in this and sub- 
sequent papers we will attempt wherever possible to "ob- 
serve" our model clusters using techniques similar to those 
employed by observers. For this paper we have chosen to 
adopt a "photometric" apjjroach. Consequently, we present 
little detailed information about the binaries in our calcu- 
lations. Binary properties will be discussed in depth in a 
"spectroscopic" companion paper (Portegies Zwart et. al. 
Paper IVb, in preparation). 

Section 2 discusses the choice of initial conditions and 
parameters for our model clusters. The results of the calcu- 
lations are presented in §3, followed in section 4, which com- 
pares our results with selected observations on a cluster-by- 
cluster basis. Section 5 our results with some previous stud- 
ies in this area. We summarize in §6. Two appendices are 
included. Appendix A gives an overview of the terminology 
used throughout the paper; Appendix B reviews the "Star- 
lab" software package and the implementation and coupling 
of its two main constituents kira (the iV-body integrator 
§B1) and SeBa (the stellar /binary evolution program§B2). 



2 INITIAL CONDITIONS 

In order to begin a simulation, a number of critical clus- 
ter parameters must be chosen: the mass, virial radius, and 
tidal radius (or its equivalent), the initial mass function, and 
the initial distributions of binary spatial density and orbital 
properties. Table 1 presents an overview of observed parame- 
ters for some star clusters with similar overall masses, stellar 
membership and half-mass radii. The ages of these clusters 
vary widely, from about 110 Myr (Pleiades) to over 1 Gyr 
(NGC3680). Their core and half mass radii suggest that 
they may be described approximately by King models (King 
1966) with dimensionless depths in the range Wo ~ 4 — 6, 
where larger Wo corresponds to higher central concentrar 
tion. 

We define our mass scale by arbitrarily adopting a 
"Hyades-Iike" model, in which the mass of the system at 
an age of 625 Myr is ^ IOOOMq. We then estimate the ini- 
tial mass of a cluster by applying a number of corrections. 
We adopt the initial mass function for the solar neighbor- 
hood described by Scalo (1986). Table 2 shows how the Scalo 
mass function evolves in time. Approximately 20% of the 
initial mass is lost by purely stellar-evolutionary processes. 
(Binary evolution complicates matters which we neglect in 
these rough numbers.) 



Table 2. Number Nx of stars of type x and total mass M, from 
the Scalo (1986) initial mass function, evolved with SeBa. Cal- 
culation was performed with a population of 100 k stars, but the 
numbers are normalized to 1 k stars. 



time [Myr] 





100 


200 


400 


600 


800 


ms 


1024 


1014.7 


1007.7 


996.7 


988.3 


982.1 


gs 





2.7 


4.3 


6.2 


6.2 


5.8 


wd 





3.0 


8.5 


17.6 


26.0 


32.7 


ns 





3.6 


3.6 


3.6 


3.6 


3.6 


M [Mq] 


624.4 


562.3 


541.5 


517.6 


501.1 


490.6 



In addition to stellar evolution, dynamical evolution of 
the star cluster and the tidal field of the Galaxy also tend 
to consume (eject) cluster stars, at a rate of about 10% per 
relaxation time (Spitzer 1987). With a current relaxation 
time for Hyades of about 400Myr (see Table 1) we estimate 
that the amount of mass lost by dynamical processes up to 
an age of 625 Myr is similar to the mass lost to stellar evo- 
lution. Adding these numbers provides a conservative lower 
limit to the amount of mass lost by the star cluster. This 
limit is conservative because we neglected the interaction 
between stellar mass loss and dynamical mass loss and the 
extra mass loss induced by interactions with giant molec- 
ular clouds. We therefore adopt an initial cluster mass of 
Mo ~ 1600 Mq. This is close to the 1800 M© for the initial 
mass of Hyades derived by Weidemann (1992). We assume a 
Scalo (1986) initial mass function, with minimum and max- 
imum masses of 0.1 M© and 100 M©, respectively, and mean 
mass (m) ~ O.6M0. Consistent with the above estimates, 
our simulations arc performed with 1024 single stars and 
1024 binaries, for a total of 3096 (3k) stars and a binary 
fraction of 50%. 

Stars and binaries within our model arc initialized as 
follows. A total of 2k single stars arc selected from the ini- 
tial mass function and placed in an equilibrium configura- 
tion in the selected density distribution (see below). We then 
randomly select half the stars and add a second companion 
star to them. The masses of the companions are randomly 
selected between 0.1 M© and the primary mass. For low- 
mass primaries, the mass ratio distribution peaJis at unity, 
whereas the distribution is flat for more massive primaries. 
Once stellar masses are chosen, other binary parameters are 
determined. Binary eccentricities are selected from a ther- 
mal distribution between and 1. Orbital separations a are 
selected with equal probability in log a with the lower limit 
set by the separation at which the primary fills its Roche 
lobe or at 1 R© , whichever is smaller. The upper limit for 
the initial semi-major axis is taken at 10^ R© (about 0.02 pc, 
Duquennoy & Mayor 1991). When a binary appears to be in 
contact at pericenter, new orbital parameters are selected. 
Table 3 gives an overview of the various distribution func- 
tions from which stars and binaries are initialized. 

We select initial density profiles from the anisotropic 
density distributions described by Heggie & Ramamani 
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Table 1. Observed and derived parameters for several open star clusters with which our simulations may be compared. Subsequent 
columns give (3) the distance to the cluster (in pc), (4) the cluster age (in Myr), (5) the half mass relaxation time (in Myr), (6) the total 
mass (in Mq), (7) the tidal radius (in pc), (8) estimate for the half mass radius (in pc), and (9) the core radius (in pc). In cases where 
the parameters (relaxation time, mass, etc.) are not accessible in the literature, we calculate it; these entries are printed in italics. In 

most cases these numbers can be calculated using Eq. A9 or Eq. AlO. Dashes and question marks indicate that wc cannot derive these 
numbers from the literature. The final two columns contain information on the cluster stellar content. The column labeled f^-.h-.t indicates 
the number of single stars, binaries and triples (separated by colons). For clusters where the numbers arc given directly by observations, 
the table gives the observed numbers of each system. If the binary fraction is derived by other methods, we give the relative fractions 
normalized to the number of single stars. The last column (ATbgg.gg.^j) gives the number of observed blue stragglers, giants and white 
dwarfs, separated by colons. -i- 



name 


ref. 


d 


t 




M 


nide 




^core 


fs:h:t 


-^bss:gs:wd 






[pc] 


[Myr] 


[Mq] 




[pc] 








NGC2516 


a 


373 


110 


220 


1000 


13 


2.9 




16:6:? 


6:4:4 


Pleiades 


b 


135 


115 


150 


~1500 


16 


2-4 


1.4 


137:60:2 


0:3:3 


NGC 2287 


c 


655 


160-200 




^ 120 


6.3 






1:0.6:? 


3:8:3 


Prafisepe 


d 


174 


400-900 


570 


1160 


12 


3.5 


2.8 


1:0.3:0.03 


5:5:11 


Hyades 


e 


46 


625 


390 


500-1000 


10.3 


3.7 


2.6 


1:0.4:0 


1:4:10 


NGC 2660 


f 


2884 


900-1200 


315 


^ 400 


9.6 


4 


1.5 


1:0.3:? 


18:39:? 


NGC 3680 


g 


735 


1450 


28 


^ 100 


4.3 


1.2 


0.6 


44:25:0 


4:17:? 



References to the literature (second column) arc: (a) Abt & Levy (1972); Daclis, J & Kabus (1989); Hawlcy ct al. (1999). (Note: we 
interpret the quoted limiting cluster radius as the half mass radius.) (b) Pinficld ct al. (1998); Raboud & Mcrmilliod, (1998); Bouvier 

et al (1998); (c) Harris et al. (1993); lanna et al (1987); Cox (1954). (d) Andricvsky (1998); Jones & Stauffcr (1991); Mermilliod & 
Mayor (1999); Mermilliod et al. (1990); Hodgkin et al. (1999). (Note: wc interpret the quoted central radius for the cluster as the half 
mass radius.) (e) Perryman et al. (1998 and references therein) Reid & Hawley (1999); (f) Prandsen et al. (1989); Hartwick & Hesser 
(1971); SandreUi et al. (1999). (g) Hawley et al. 1999 ; Nordstrom et al. (1997); Nordstrom et al. (1996), Data on numbers of white 
dwars was taken from Anthony- Twarog (1984) for Praesepe, from Koester & Reimers (1987) for NGC 2287 and from von Hippel (1998) 

for the other clusters. 



Table 3. Initial conditions for the stellar and binary popula- 
tion. The first column gives the parameter, the second and third 
columns give the symbol and the distribution function, followed 
by the lower and upper limits adopted. 



parameter 
primary mass 
secondary mass 
orbital separation 
eccentricity 



M 
m 

a 

e 



function 

PiM) = Scalo (1986) 
P(m) = constant 
P(a) = 1/a 
P{e) = 2e 



limits 
lower upper 
0.1 Mq 100 M0 
0.1 Mq M 
RLOF 10^ Rq 
1 



(1995) with Wo = 4 and Wo = 6, and refer to these 
models as W4 and W6, respectively throughout this pa- 
per. The Hcggic-Ramamani models arc derived from King 
(1966) models, but take into account the velocity anisotropy 
and non-spherical shape of the critical zero- velocity (Jacobi) 
surface of the cluster in the Galactic tidal field. (The clas- 
sical King models have spherical boundaries.) Within the 
half mass radius, the Heggie-Ramamani models are quite 
isotropic. 

All models are started with a virial radius of rvir = 2.5 
pc. For our adopted parameters, the initial cluster dynamical 
time scale is then thm = (GM/rvir'^)"^^^ ~ 1.5 Myr. Each 

cluster is assumed to precisely fill its Jacobi surface at birth. 
(Expressed less precisely, we could say that the limiting ra- 



dius of the initial King model is equal to the Roche radius 
of the cluster in the Galactic tidal field.) Given the Oort 

constants in the solar neighborhood, we find that the mod- 
els with Wo = 6 are somewhat farther ('--^12.1 kpc) from the 
Galactic center than is the Sun, while a model with Wo = 4 
Is slightly closer (~6.3kpc). 

For a total cluster mass of 1600 Mq , the Lagrange points 
of our two standard clusters lie, respectively, at 14.5 pc 
(Wo = 4) and 21.6 pc (Wo = 6) from the cluster center. 
A star is removed from a simulation when its distance from 
the cluster's density center exceeds twice the distance from 
the center to the first Lagrangian point. 

Table 4 reviews the adopted parameters and initial con- 
ditions of our models. In order to improve statistics, we per- 
formed four calculations (labeled I through IV) for each set 
of initial conditions. A fifth run is performed as a pioneering 
study for each set of initial conditions but they are termi- 
nated at 400 Myears. 



3 RESULTS 

We now discuss the "photometric" properties of our model 
clusters. As mentioned above, we defer the discussion of 

"spectroscopic" properties, including details on the various 
types of binaries found in our simulations, to Paper IVb. 
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Table 4. Initial conditions and parameters for the selected models. The columns give the model name, initial mass (M©), number of 
stars, dimensionless central depth of the potential well (Wo), the distance from the Galactic center (kpc), the initial relaxation and half 
mass crossing times (both in Myr) , , ry, and , the distances from the cluster center to the Ja<;obi surfa<;e (so rx is simply the Jacobi 
radius, rj; see Appendix A), and the virial, half mass and core radius (all in parsec). 



name 


M 


N 


Wo 


^Gal 






''Jacob 








''core 




[Mq] 






[kpc] 


[Myr] 






[pc] 




[pc] 




W4 


1600 


3k 


4 


6.3 


109 4.07 


14.5 


9.7 


7.2 


2.5 


2.14 


0.83 


W6 


1600 


3k 


6 


12.1 


102 4.15 


21.6 


14.4 


10.8 


2.5 


2.00 


0.59 



Table 5. Overflew of the variation in parameters for model W6- 
III. The realization of the initial mass function was identical for 
all these models. Other parameters are as in Tab. 4. For each of 
these calculations some feature of starlab was switched on (+) 
or off (-); kira(the W-body integrator), SeBa (the stellar and bi- 
nary evolution model) and wether the calcualtion started with 
primordial binaries. 





Starlab 
kira SeBa 


Primordial 
binaries 


upper solid 


+ 


+ 


lower solid 


+ 


+ 


dashes 


+ + 


+ 


dash-dots 


+ + 





3.1 Global properties 

Figure 1(a) shows the mass of the cluster as a function of 
time for several models. In this figure, the "total mass" of 
a model is simply taken to be the sum of the masses of 
all stars remaining within the Af-body system. This over- 
estimates both the bound mass of the system and proba- 
bly also the mass that would be derived by observers, and 
therefore provides a firm upper limit to the "true" mass. 
The dashed line gives the results from model W6-III. The 
two dotted lines show the mass evolution of two of the W4 
models (upper line: model W4-III, lower line: W4-IV), il- 
lustrating the run-to-run variations in dynamical evolution 
(which are mainly the result of an initial offset between the 
masses of the two models, caused by different random seeds). 

Table 5 indicates the differences between the various line 
styles in Fig. 1(a). All lines listed are computed using exaclty 
the same stellar masses as in model W6-III. 

The two solid lines in Figure 1(a) show (upper line) the 
total mass of an A^-body system without stellar evolution, 
but otherwise with the same initial conditions as run W6-III, 
and (lower line) the total mass of stars, excluding stellar dy- 
namics but including mass loss by stellar evolution (with the 
same initial stellar masses as in model W6-ni; see Table. 5 
for an review of the line styles) . Mass loss in the absence of 
stellar evolution is not linear with time, as one expects from 
a equal mass AT-body system; the presence of a mass func- 
tion causes heavier stars generally to be lost later, resulting 
in larger mass loss at later time. However, the loss rate of 



stars is roughly constant, and the curvature of the upper 
solid line in Figure 1(a) demonstrates the strong effects of 
mass segregation and cluster dynamics. 

The actual variation of the cluster mass (under the com- 
bined effects of stellar mass loss and dynamical evolution) is 
larger than the sum of the two separate effects by about a 
factor two. For the selected initial conditions, the interplay 
between stellar evolution and stellar dynamics is especially 
important during the later stages of the evolution. The pres- 
ence of primordial binaries has little effect on the evapora- 
tion rate of the clusters. The dash-dotted line gives the total 
mass of the same model, but without binaries — all stars in 
the initial mass function (including binary secondaries) are 
redistributed with the same density profile as model W6-III. 
As noted by McMillan & Hut (1994), primordial binaries 
have little effect on the overall rate at which mass is lost 
from the cluster. 

Figure 1(b) shows the masses of two model clusters 
(W6-HI and W4-H), with different criteria for the limiting 
radius of the stellar system. The solid lines give the total 
mass in the Af-body system (see also Figure l[a]), the dot- 
ted lines the total mass within the Jacobi surface, and the 
dashed lines the total mass within the Jacobi radius of the 
cluster center, as seen by an observer looking along the y- 
axis. The total mass in stars in the Af-body system overesti- 
mates the cluster's mass; the mass within the zero velocity 
surface may give a better measure of the mass one would ob- 
serve in a real situation. The cluster mass within the Jacobi 
radius, viewed along the x-axis, always lies between the solid 
and dotted lines; the corresponding mass viewed along the 
z axis lies between the dotted and the dashed lines. These 
trends are found in all models studied. 

Figure 2(a) presents the observational equivalent of Fig- 
ure 1(b), showing the total Adv magnitude of the W6 models 
at various times. The spread in My is caused by the intrin- 
sic differences bctwccu runs — initial conditions, run-to-run 
variations and fluctuations due to the small numbers of gi- 
ants, which dominate the total magnitude. Note that, for 
technical reasons, the output intervals are not the same for 
the four calculations shown here. 

Figure 2(b) plots the integrated B — V color of the 
cluster as a function of time. The initial color of all mod- 
els is confined to a small range between B — V ^ —0.15 
and —0.28. but rapidly grows to larger values with a larger 
spread: B - V ~ 0.1-0.4 at 50 Myr and B - F ~ 0.5-0.8 
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Figure 1. (a) Mass (in Mq) as a function of time (in Myr) for various models. The solid lines show the time-dependence of the total 
mass with (upper curve) stellar evolution suppressed, and (lower curve) without dynamical evolution. The dashed line is the total mass 
of model W6-III. The dot-dashed line is the total mass of a calculations with the identical mass function as in model W6-III, but all 
stars were single and redistributed in an identical density distribution. This line therefore gives the difference in a calculation with or 

without primordial binaries, but with all the other effects takes the same. The dotted lines represent the total mass of models W4-II 
(upper) and W4-IV(lower). (b) Various definitions of the cluster mass (in M0) as functions of time for model W6-III (upper set of 3 
lines) and model W4-II (lower set of lines). Solid lines represent the total mass in the Af-body system, dotted lines the mass within the 
Jacobi surface, and dashed lines the total mass seen in projection within a distance rj of the cluster center. 



+ s ♦ t„ . 
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Figure 2. (a) Integrated My magnitude as a function of time of the four W6 models (plus, circle, square, and diamond give the data 
for the models I-IV) . The variation between the points at similar times within a single calculation is comparable to the visible spread in 
the data, (b) Integrated B — V color as a function of time for the W6 models. 



at 500 Myr. As in Figure 2(a), the intrinsic spread is caused 
by the presence of a relatively small number of giant stars. 
The increase in the color index indicates that the cluster 
gets redder with age, which is mainly duo to the loss of the 
massive blue stars and the formation of red giants. The color 
variation of the W4 models is similar. 

Figures 3 shows how the radii of models W4 and W6 
evolve with time. All stars in the iV-body system are taken 
into account in calculating the Lagrangian radii. The outer 
radii therefore expand much more than they would do if 
only stars within the Jacobi surface were considered. Note 
the absence of any discernible core collapse in either case. 
There is a slight, barely noticeable, core contraction between 



t = 100 and 150 Myr for the W6 models, and somewhat ear- 
lier for the W4 models, but neither is very deep. This shallow 
core contraction phase demonstrates the importance of stel- 
lar mass loss and, to a lesser extent, binary heating to the 
dynamical evolution of these systems; mass segregation for 
example, which also plays an important role here (see Fig. 5). 
A comparable model with primordial binaries and without 
stellar evolution would experience core collapse (McMillan 
et al. 1990, 1991), even in the presence of a Galactic tidal 
field (McMillan & Hut 1994). 

In contrast to the W6 radii shown in Figure 3(a), the 
Lagrangian radii of model W4-II (Figure 3[b]) expand for 
about 300 Myr, and subsequently shrink. The decrease in 
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Figure 3. Lagrangian radii (in pc) as functions of time. The data in (a) represent the mean of the four W6 runs, those in (b) model 
W4-II. Prom top to bottom, the radii contain the following mass fractions: 75% (dots), 50% (upper solid), 25% (dashes) and 5% (lower 
solid). 



Lagrangian radii at late times is an indicator of cluster evap- 
oration. As the cluster dissolves in the Galactic tidal field, 
its tidal radius decreases, accelerating the dissolution and 
causing the Lagrangian radii to decrease. The W6 clusters 
show the same behavior, but at somewhat later times. Wo 
show the result of a single W4 model to illustrate the intrin- 
sic fluctuations within a single AT-body run. 

Finally, Figure 4 shows the half-mass relaxation time 
(Eq. A9) as a function of time for models W6-III (solid line) 
and W4-II (dashed line) . Note that the relaxation time peaks 
around the cluster's "half-life" epoch — 750 and 400 Myr for 
models W6-III and W4-II, respectively. Consequently, esti- 
mates of the present-day relaxation time of observed open 
clusters may provide poor indicators of the dynamical age 
of the stellar system. The often-used argument that a star 
cluster is barely older than its relaxation time and therefore 
cannot be dynamically evolved is clearly in error for the ma- 
jority of star clusters (see also McMillan & Hut 1994). 

3.2 Mass segregation 

The effect of mass segregation is clearly visible in Figure 
5(a), which shows the mean mass (m) of stars within the 5%, 
25%, 50% and 75% Lagrangian radii as functions of time, 
averaged over the four W6 models. The initial increase in 
the mean mass within the inner 5% Lagrangian radius is the 
result of mass segregation. The value of (m) in the cluster 
center decreases again after about ICQ Myr, when the most 
massive stars leave the mean sequence and lose most of their 
mass on the Asymptotic Giant Branch. For the remainder 
of the calculation {m) stays more or less constant in each 
Lagrangian zone, but with a significantly higher value in the 
inner zones. 

Figure 5(b) shows the evolution of the mean mass (m), 
in model W4-II. Mass segregation in this model proceeds on 
a longer time scale than in the W6 models, but the cluster 
dissolves on a shorter timescale. Near the end of the cluster 




500 750 1000 

t [Myr] 

Figure 4. Half-mass relaxation time as a function of time for the 
models W6-III (solid) and for model W-HI (dashes). The error 
bars indicate the computed relaxation times for the observed star 
clusters from Table 1. (Confidence intervals are not listed in the 
table.) 



lifetime the mean mass in the outer regions increase, caused 
by the preferential loss of the lower mass stars by evapora- 
tion. The more massive stars have greater difficulty to climb 
out of the potential well. 

Figure 6 shows the mean mass to light ratio for the W6 
models. Except for the first few million years, the mass-to 
light ratio in the cluster core is always substantially smaller 
(and noisier) than that in the halo, as mass segregation 
causes the most massive (i.e. brightest) stars to sink rapidly 
to the cluster center. The general increase in mass-to light 
ratio with time is the result of stellar evolution, the loss of 
lower- mass stars by tidaJ stripping compensates somewhat. 
The occasional dips in the mass-to-light ratio are caused by 
individual red giant stars; most of these dips occur in the 
core, where mass segregation causes the giants to accumu- 
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Figure 5. (a and b) The mean mass (m) of model W6-III (a) and W4-II (b), as functions of time. The data have been smoothed over 
time intervals of 8.8 Myr. Prom top to bottom, the lines represent the mean mass within the 5% (solid), 25% and 50% (dashes), and 75% 
(dots) Lagrangian radii (see Fig. 3 for the corresponding radii.) 
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Figure 6. The mean mass to light ratio for models W6-I to IV, 
within the 5% (solid) and 75% (dots) Lagrangian radii. 



late. We can reduce this "noise" considerably by averaging 
over the four W6 calculations, but the cflFcct remains visible. 

Figure 7(a) shows the radial stellar distribution in the 
W6 models at various epochs. The cluster expands as it ages. 
The binaries (dotted lines) closely follow the distribution 
of the single stars for the first 100 Myr, but become more 
centrally concentrated at later epochs. Mass segregation is 
also clearly visible if we compare the radial distributions of 
low mass (faint) stars with the more massive (bright) stars 
(Figure 7[b]). Stars with L > 0.5 Lq are clearly more cen- 
trally concentrated than the mean cluster star, while giants 
(although there arc only a few) are even more strongly con- 
centrated in the cluster center. 

Mass segregation can also be observed in the cluster's 
mass and luminosity functions. Figure 8(a) shows global 
mass functions for all single stars and binary primaries for 
the W6 models at birth and at t = 600 Myr. The global mass 
function of the cluster is affected only slightly by stellar evo- 



lution and mass segregation. liowevcr, the mass function at 
t = 600 Myr for stars in the inner part of the cluster (dot- 
dashed line) is clearly different from the global mass function 
at that time. 

The white dwarfs arc more centrally concentrated than 
stars with luminosity > 0.5 L© and slighly less concentrated 
that the giants (see Fig. 7[b]). This is caused by the progen- 
itors of the white dwarfs, the giants, being centrally concen- 
trated while after their envelepes are shed they have masses 
comparable to the the mean cluster stars. Segregating out- 
wards takes more time than sinking inwards and at the same 
time more white dwarfs are produced in the cluster center 
(sec the end of this section for more detailes). 

Figure 8(b) shows the luminosity functions (in AIv) for 
stars (including binaries) in models W6 at zero ago and at 
600 Myr. Note that the luminosity function at the later time 
shows a larger faction of bright stars. The reason is two fold: 
1) stellar evolution has removed only the most massive stars 
by this time, while mass segregation has concentrated the 
remaining massive (bright) stars in the cluster core, at the 
same time causing the lower mass stars to escape and 2) 
heaviest stars turn into giants, which for older (lighter) stars 
arc much brighter than main-sequence stars, whereas for 
younger (heavier) stars, the difference in brightness between 
giant and main-sequence star is much smaller. 

The W4 models more strongly affected by mass seg- 
regation (not shows but see sect 4). In part, this is caused 
by the more rapid evaporation of these models compared 
to the W6 models. This is consistent with the findings of 
Takahashi & Portegies Zwart (2000), who noted that clus- 
ters which are close to complete disruption contain a higher 
fraction of high-mass stars. 

3.3 Hertzsprung— Russel diagrams 

Figure 9 shows a time sequence of Hertzsprung-Russel dia- 
grams for model W6-111. The youngest Hertzsprung-Russel 
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Figure 7. (a) Cumulative distribution of single stars (solid lines) and binaries (dotted lines), averaged over the four W6 models, at 
(upper left to lower right) t = 0, 100, 200, 600, and 800 Myr. (b) Cumulative distribution of various stellar populations for the four W6 
runs at an age of 600 Myr. Solid line: the distribution of all stars (see also the lower solid line in Figure 7[a]); dash-dotted line: stars 
with L > 0.5 Lq; dashed line: white dwarfs; dotted line: (sub)giants. 
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(a) (b) '''' 

Figure 8. mass function (a) and luminosity function (b) for model W6. The dotted lines give the f = distributions and the dashed 
lines give the distributions at t = 600 Myr (averaged over all models). The dash-dotted lines are the mass (luminosity) functions for 
stars within 3.4 pc (the half-mass radius for W6 models at t = 600 Myr) of the density center (for [b] in projection, as viewed along the 
X axis). The solid line gives the mass fucntion of the W6 models within 3.4 pc (projected again for [b]) at an age of llOOMyear. 



diagram (200 Myr) already shows a white dwarf sequence. 
Note the densely populated "binary sequence" ~0.75 magni- 
tudes above the zero-age main-sequence. One of the objects 
in the middle panel (close to but just above the turnoff) is 
a blue straggler; the other two arc binaries (sec also Figure 
11). The objects immediately to the left of the main sequence 
(the two points in the 600 Myr diagram at _B — V ~ 1.18 
and the single point at B — V ^ 0.5, V = 8) are binaries 
containing a mass-transfer remnant (a helium star) and a 
main-sequence star which has accreted part of its compan- 
ion's envelope. Farther to the blue (between B — V = 
and 0.8), but to the right of the white dwarf sequence, are 
binaries consisting of a white dwarf and a low-mass main- 
sequence star (several are seen in the 600 Myr and 1100 Myr 
diagrams) . In the bottom panel a break and discontirmity in 
the zero-age main-sequence is visible near B — V = 0.4 and 



at V ~ 4. This is an artifact of the stellar evolution fitting 
formulae given by Eggleton, Fitchet &Tout (1989) and ap- 
pears when the envelope of a main-sequence star becomes 
convective. 

Figure 10 shows Hortzsprung-Russel diagrams for the 
irmcr, middle and outer regions of the combined W6 mod- 
els at an age of about 600 Myr, and illustrates the effect 
of mass segregation. Each diagram contains about 2000 ob- 
jects. The slight "fuzziness" near the main-sequence turnoff 
is the result of variations in output times between individ- 
ual simulations, which cause the combined diagram to have 
a small spread in stellar ages. About one quarter (one run) 
of the stars come from a slightly younger cluster with an age 
of about 550 Myr. 

The Hertzsprung-Russel diagrams of the inner and 
outer parts of the cluster show significant difi'erences due 
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Figure 9. Hertzsprung-Russel diagram of all stars in model W6-III at ages of 200 Myr U — B in the left panels and _B — y in the right 
panels (upper panels; 1983 objects), 600 Myr (middle; 1603 objects), and 1100 Myr (lower panels; 1009 objects). 



to mass segregation. Most are a consequence of the evolv- 
ing binary population, and will be discussed in more de- 
tail in Paper IVb. The inner HRD has a clear excess of 
(sub)giants and white dwarfs relative to the HRD at the 
half mass radius or that in the halo. Also, the turnoff region 
is more heavily populated in the inner HRD than in the 
others. Striking also is the lack of a clear binary sequence in 



the outer Hertzsprung-Russel diagram. The bottom of the 
main-sequence is less clearly affected by mass segregation. 



3.4 Blue Stragglers 

The (small) numbers of blue stragglers do not depend 
strongly on the particular region of the cluster under study. 
We count four blue stragglers in the inner Hertzsprung- 
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Figure 10. Hertzsprung-Russel diagrams of the combined W6 
models at an age of about 600 Myr. The upper panel shows the 
innermost (non-projcctcd) 4 pc (2004 objects), the middle panel 
stars between 4 and 9 pc from the cluster center (2658 objects), 
and the bottom panel stars more than 9 pc from the center (2039 
objects). 



Russel diagram, and one and two in the middle and outer 
frames of Figure 10, respectively. These numbers are fairly 
typical of our simulated clusters and also quite typical for 
the numbers observed (sec Tab. 1). Fig 11 presents a graphi- 
cal representation of the blue stragglers in model W6-III. (A 
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Figure 11. Blue stragglers in model W6-III. The solid curve 
gives the turnofl mass (in M©) as a function of time (in Myr). 
The horizontal lines represent the tracks of the blue stragglers in 
model W6-III. The tracks start when the star is rejuvenated (see 
Appendix B), and stop when the blue straggler leaves the main 
sequence. 



main-sequence star is identified as a blue straggler as soon 
as its mass exceeds the turnoff mass for that epoch.) 

Most blue stragglers are the result of mass transfer 
in a close binary. In about half of all cases (37 out of 76 
blue stragglers formed in all calculations performed), the 
mass transfer is unstable, leading to a merger. Blue strag- 
glers formed from a stable phase of mass transfer are gen- 
erally accompanied by a white dwarf or helium star (the 
young remnant of mass transfer, see Figure 11), causing the 
blue stragglers to lie slightly blucward of the turnoff in the 
Hertzsprung-Russel diagram (see §3.3). 

Three blue stragglers formed via collisions in which a 
third star interacted with and became bound to binary, lead- 
ing to a collision between the binary components. The orbit 
of such a blue-straggler binary is generally quite eccentric, 
and the companion to the blue straggler is most likely to 
be a main sequence star. Observing a blue straggler in an 
elliptical orbit around a stellar (main-sequence) companion 
would provide strong evidence for such dynamical interac- 
tions in star clusters (see Portegies Zwart 1996). 

In none of our calculations a blue straggler with more 
than twice the turn off mass was formed i.e., there where 
no collisions between throe or more stars. A discovery of a 
blue straggler with a mass more than twice the turn-off mass 
would provide strong evidence for effects of stellar dynam- 
ics, though one could imagine a primordial triple to get into 
a common-envelope situation in which all three stars spi- 
ral in to a triple merger. Portegies Zwart ct al. (1999) find 
runaway collisions between more than two stars in their sim- 
ulations of dense stars clusters without primordial binaries. 
Thier results, however, are applicable for a different range 
of initial conditions, as they studied the dence and young 
central star clusters R 136 in the 30 Doradus region of the 
Large Machelanic Could. 
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Many blue stragglers experience mass transfer or a col- 
lision long before actually being classified as blue stragglers 
by our criterion (i.e. exceeding the turnoff mass). In most of 
these cases, one or more phases of mass transfer (stable or 
unstable, or even accretion from the stellar wind of a com- 
panion) has rejuvenated one of the stars in a close binary 
system (see Appendix B). As the cluster ages, the star re- 
mains behind on the main sequence, and eventually becomes 
identifiable as a blue straggler (see also paper II) . This is il- 
lustrated in Figure 11 (the two long tracks with m ~ 2 and 
the track near m 3) . 

A blue straggler which was rejuvenated long ago may 
show no trace of the event that caused its rejuvenation. 
Apart from residing above the turnoff, the star may appear 
completely normal; anomalous atmospheric abundances will 
have had sufficient time to mix with the stellar interior. In 
addition, if the blue straggler is rejuvenated only a little, 
the maximum distance on the Hertzsprung-Russol diagram 
between the cluster turnofi' and the blue straggler will be 
very small; such a star may remain unidentified as a blue 
straggler. This may happen if mass transfer is unstable but 
does not lead to a merger, or if a binary is too wide for 
Roche-lobe overflow, and the blue straggler is rejuvenated 
by accreting a small portion if its companion's stellar wind. 

The lifetime of a blue straggler depends on the epoch 
at which it formed. Blue stragglers that form later are gen- 
erally products of lower-mass stars, and tend to live longer 
than blue stragglers that formed early in the evolution of 
the stellar system. 
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3.5 Isophotes 

Figure 12 shows a series of isophotes, as seen from vari- 
ous directions, for model W6-III at an age of 600 Myr. The 
Galactic center is located to the —x direction (at a distance 
of about 12.1 kpc — see Table 4), and z points toward the 
Galactic north pole. While the cluster is barely flattened at 
birth, § by 600 Myr the cluster is significantly flattened by 
the Galactic tidal field. As expected, the flattening is great- 
est along the z axis, and also noticeable in the y direction. 

Figure 13 shows images of model W6-III at three differ- 
ent moments in time. The images are created using a ray- 
tracing technique. 

3.6 Escaping stars 

Stars escaping from the cluster arc lost primarily near the 
first and second Lagrange points. Figures 14 and 15 show the 
positions and projected velocities of the first 100 escapers 
from the system, and the 100 stars which escaped between 

§ This is simply a consequence of the fact that only the outer- 
most parts of the cluster, near the Jacobi surface, show significant 
flattening, and these are initially very sparsely populated. Only 
when cluster evolution drives many stars out to the Jacobi radius 
does the flattening become readily apparent. 
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Figure 12. Isophotes in My for model W6-III at a cluster age of 
600 Myr. The three panels present views along the three coordi- 
nate axes. The 10 (sub)giants are plotted as dots with size propor- 
tional to magnitude (and are excluded from the isophotes). The 
brightest region in the continuum plot has a surface brightness 
of -2.2 mag/pc^^. Contours are plotted at -0.86 magpc^^, 0.34 
magpc"^, 1.6mag pc^^, 2.8 magpc"^, 4.1 magpc"-^, and 7.8 
magpc"^. Stars are assigned a Gaussian point spread function 
with a dispersion of 0.35 pc. 
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Figure 13. Visualization of model W6-III at zero age (top im- 
age), at an age of 622 Myr and at an age of 1512 Myr. Images 
were created using a ray-tracing technique. 

t = 550 and t — 650 Myr. The left and right panels show, 
respectively, projections onto the x — z and x — y planes. 
The first and second Lagrange points lie on the x-axis, at 
distances of ~20 {t = 0) and ~17 (i 600) pc from the 
cluster center. The small overall rotations of escapers evident 
in the X — y projections arc consequences of the Coriolis 
force acting on stars in the rotating frame of reference in 
which we perform the simulations. The high-speed escapers 
with roughly isotropic velocities in Figure 14 are escaping 
neutron stars, which receive high kick velocities on their 
formation. They arc absent in Figure 15, as the cluster is 
by that time too old for supernovae to occur (except for 
type la supernovae). 

The main differences between Figs. 14 and 15 are (1) 
the considerably larger spread in velocities, (2) the larger 
extent in z of the region over which stars arc lost, and (3) the 
higher speeds of escaping stars at the earlier epoch. These 
differences are readily explained by a combination of effects; 
the evolution of the cluster in the Galactic tidal field, the 
presence of primordial binaries and the formation of neu- 
tron stars. As the cluster ages it becomes less massive and 
the Galaxy's gravitational pull becomes relatively stronger. 
The tidal radius shrinks and the cluster velocity dispersion 
decreases, so the speed of escaping stars and their distances 
above or below the Galactic plane also decrease. The older 
cluster also lacks massive stars and no stars are ejected via 
supernova explosions. The shallow core collapse during the 
first 100 Myr results in increased binary activity, which also 
contributes to the higher stellar ejection speeds at the earlier 
time. 

3.7 Stellar populations 

Tables 6 and 7 present, for several cluster ages, the num- 
bers of single stars and binaries by generic stellar types. 
The numbers are averaged over the calculations performed 
for each set of initial conditions, chosen with a different ran- 
dom seed from the same probability distribution. Table 8 
gives the same data for a population of evolving binaries 
without dynamics, calculated using SeBa (see Appendix B). 

Overall, the evolutionary differences in the populations 
of single stars and binaries between models W6 and W4 are 
quite small. Clearly, as already noted, the W4 clusters evap- 
orate more rapidly (Figure 1), resulting in a generally more 
rapid decrease in the numbers of both stars and binaries. 
A more significant difference between Tables 6 and 7 is the 
larger rmmbers of white-dwarf binaries in the W4 models 
compared to other stellar types. Table 8 presents the stellar 



Table 6. Stellar and binary types in the W6 models at various 
times. Binaries which contain two main-sequence stars are iden- 
tified as (ms, ms), (ms, gs) contain a main-sequence star and a 
giant, (gs, gs) contains two giants, (gs, wd) contains a giant and a 
white dwarfs and (wd, wd) contans two white dwarfs. A bracket 
indicates the binary component which fills its Roche-lobe and is 
in a state of mass transfer to its companion star, binaries with 
neutron stars or black holes are omitted. The bottom row gives 
the binary fraction. 



# runs: 


5 


5 


5 


5 


4 


4 


time [Myr]: 





100 


200 


400 


600 


800 


ms 


1024 


939 


915 


830.4 


914.5 


508.5 


gs 





3.2 


3.8 


5.0 


6.0 


4.5 


wd 





4.8 


9.4 


21.2 


33.3 


26.3 


(ms, ms) 


1024 


688.6 


644.2 


603.0 


520.5 


438.0 


[ms, ms) 





4.2 


3.4 


3.2 


3.3 


2.5 


(ms, wd) 





1.0 


2.4 


5.4 


8.8 


9.5 


(gs, ms) 





0.4 


2.4 


1.6 


4.3 


2.8 


(gs, gs) 





0.0 


0.4 


0.2 


0.3 


0.3 


(gs, wd) 





0.4 


0.4 


0.6 


1.5 


2.0 


(wd, wd) 





0.0 


1.2 


3.6 


6.0 


7.5 


/bin 


0.5 


0.42 


0.41 


0.42 


0.36 


0.46 



Table 7. Stellar and binary types in the W4 models. 



# runs: 
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5 
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5 
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2 


liuK! iAl\"r^: 





100 


200 


100 


(iOO 


SOD 


ms 


1024 


998.8 


938.8 


771.2 


604.5 


164.0 


gs 





2.8 


4.2 


7.4 


9.8 


7.5 


wd 





4.2 


9.4 


20.6 


33.5 


31.0 


(ms, ms) 


1024 


502.8 


600.6 


453.0 


333.0 


187.5 


[ms, ms) 





4.8 


4.1 


2.8 


2.3 


1.1 


(ms, wd) 





1.8 


2.2 


2.4 


5.5 


6.0 


(gs, ms) 





0.6 


2.4 


1.2 


2.5 


3.0 


(gs, gs) 





0.0 


0.0 


0.2 


0.0 


1.0 


(gs, wd) 





0.4 


0.6 


0.4 


1.0 


0.0 


(wd, wd) 





0.2 


1.6 


3.0 


4.0 


5.0 


/bin 


0.5 


0.34 


0.39 


0.37 


0.35 


0.50 



and binary properties of an evolving population of isolated 
binaries. The differences between these binaries and the dy- 
namically evolving population is considerable. Comparing 
Table 8 with Tables 6 and 7 reveals that the dynamically 
evolving populations are enhanced in both giants and white 
dwarfs; the effect is stronger in the W4 models because of 
the enhanced escape of the lighter stars. 

Table 9 gives the fractions of various types of stars and 
binaries in the dynamical calculations, relative to the cor- 
responding numbers from the population synthesis studies. 
The latter are normalized to the same numbers of single 
main-sequence stars and main-sequence binaries as in the 
initial dynamical calculations. (This normalization is cm- 
ployed here to show trends which are hard to see in the 
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Figure 14. Vector diagram of the first 100 (t < 218 Myr) stars escaping from model W6-III. Projections onto the x — z plane and the 
X — y plane are shown. A velocity scale is shown in the upper right corner. 
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Figure 15. As for Figure 14, but for the 100 stars escaping from model W6-III between t = 550 Myr and t = 650 Myr. 



Tables 6, 7 and 8.) Neutron stars and black holes are omit- 
ted from the comparison because the differences are directly 

obvious: Neutron stars escape from open star clusters, but 
they are retained in the non-dynamical models. Black holes 
are omitted because of their small numbers. 

Although the numbers of giants and white dwarfs are 
also small, a trend is clearly visible: Single white dwarfs and 
dwarfs in binaries are overrepresented at later stages in the 
dynamical calculations, especially for the W4 models. The 
basic reason for this overabundance of white dwarfs and gi- 
ants is their larger mass, which makes them more likely to 
be retained by the cluster. White dwarfs have a very compli- 
cated evolution within the stellar system, their progenitors 
being among the most massive objects while on the main 
sequence and the giant branch, but the white dwarfs hav- 
ing masses comparable to the mean once their envelopes axe 
lost. White dwarfs axe therefore preferentially formed in the 
cores of star clusters. Once the white dwarf is formed, it is 
hard to extract it from the core. These effects are more pro- 
nounced in smaller clusters (trix ^ IGyr). The W4 models 



retain more white dwarfs than the W6 models because the 
latter relax on a longer time scale and evolve dynamically 
less rapidly than the former. 



4 COMPARISON WITH OBSERVATIONS 
4.1 The Pleiades 

Figure 16 shows the Af/-magnitude luminosity function for 
the inner part of the Pleiades cluster (Hambly & Jameson 
1991) and compares it with our model luminosity functions 
at 100 Myr. 

The best fit between the observed and model luminosity 
functions is obtained for the stars within the half mass ra- 
dius. This suggests that some mass segregation has already 
occured in this clusters. Raboud & Mermilliod (1998) also 
find evidence for mass segregation in this cluster. Our lu- 
minosity function has too many bright stars and to make a 
reasonable fit we have to exclude stars with Afj < 4.5 from 
the sample. We are not sure why this is the case, but argue 
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Table 8. Stellar types from population synthesis studies of 
1.5 X 10® binajries. The numbers of binaries is renormalized to 
1024 because this is the number of primordial binaries in e3x:h. of 
our dynamical calculations. The evolution of population of 1024 

single stars was presented in Tab. 2. Note that the dynamical 
models were performed with 1024 primordial binaries and 1024 
single stars. We added the extra class of binaries which includes 
a neutron star or a black hole, such binaries are omitted in the 
dynamical models due their small number. 



time [Myr]: 



ms 





0.79 


0.42 


0.92 


0.34 


gs 





0.44 


0.51 


0.73 


1.10 


wd 





0.87 


1.96 


3.57 


4.99 


ns/bh 





4.37 


4.42 


4.42 


4.43 


(ms, ms) 


1024 


985.14 


976.72 


965.18 


956.62 


(ms, gs) 





1.81 


3.01 


4.45 


4.82 


(ms, wd) 





1.89 


4.35 


8.48 


12.25 


(ms, ns/bh) 





0.10 


0.08 


0.05 


0.05 


(gs, gs) 





0.13 


0.25 


0.24 


0.27 


(gs, wd) 





0.27 


0.78 


1.33 


1.59 


(gs, ns/bh) 





0.01 


0.01 


0.01 


0.00 


(wd, wd) 





0.52 


1.98 


4.55 


6.78 


(wd, ns/bh) 





0.27 


0.25 


0.23 


0.22 


(ns/bh, ns/bh) 





0.11 


0.10 


0.10 


0.09 



Table 9. Relative numbers of stars and binaries in the dynamical 
models, as fractions of the numbers found in the non-dynamical 
population synthesis studies. The normalization is such that the 

dynamical and non-dynamical calculations contain equal numbers 
of single main-sequence stars and main-sequence binaries. Num- 
bers greater than 1 indicate excesses of those stellar type in the 
dynamical calculation; numbers less than 1 represent depletion. 



Normalized data for the W6 models 
time [Myr]: 100 200 400 600 800 



ms 


1 


1 


1 


1 


1 


1 


(ms, ms) 


1 


1 


1 


1 


1 


1 


gs 





1.0 


0.9 


0.9 


0.9 


1.3 


wd 





1.4 


1.0 


1.3 


1.2 


1.3 


(ms, gs) 





0.3 


1.1 


0.6 


1.2 


1.3 


(ms, wd) 





0.8 


0.8 


1.0 


1.0 


1.3 


(gs, gs) 





0.0 


2.3 


1.3 


1.3 


2.4 


(gs, wd) 





2.1 


0.7 


0.7 


1.3 


2.6 


(wd. \v<[) 





D.O 


0.<) 


1.2 


1.2 


1.8 




Normalized data for the W4 models 


time [Myr]: 





100 


200 


400 


600 


800 


ms 


1 


1 


1 


1 


1 


1 


(ms, ms) 


1 


1 


1 


1 


1 


1 


gs 





1.0 


1.0 


1.5 


2.3 


6.3 


wd 





1.2 


1.0 


1.3 


1.9 


4.6 


(ms, gs) 





0.6 


1.2 


0.5 


1.4 


3.3 


(ms, wd) 





1.7 


0.8 


0.6 


1.2 


1.8 


(gs, gs) 





0.0 


0.0 


1.6 


0.0 


2.2 


(gs, wd) 





2.6 


0.1 


0.6 


1.7 


0.0 


(wd, wd) 





0.9 


0.5 


1.3 


1.6 


2.7 



0.31 
1.11 
6.55 




4.43 

949 igure 16. Cumulative luminosity function in M/-magnitudes. 
4 4yhe open circles (o) with (Poissonian) error bars show the ap- 
gparent luminosity function for the Pleiades (Hambly & Jameson 
Q q4991). The corrected absolute luminosity function, assuming dis- 

Q 2lanec modulus m—M = 5.5 (Gatewood et al. 1990) is indicated by 
X.g§lled circles (•). Both luminosity functions are corrected for the 
0(1)50 stars brigtcr than nij = 10.5 (Pinfield ct al 1998). The dotted 
g g^ne is the initial luminosity function for all stars in the models (as- 
0.2iuming that binaries are unresolved). The solid and dashed lines 
Q Qghows the liuminosity function for all models at t = 100 Myear 
within a projected (onto the y-z plane) 25% Lagrangian and the 
half mass radius of the cluster. The dash-3dotted line give the 
luminosity function for the cluster stars in the outer 90% La- 
grangian radius. 



that the brightest stars were possibly overexposed in the ob- 
servations, and may therefore have been omitted from the 
observed luminosity function. 

The Pleiades is flattened, with an observed ellipticity 
e = {1 — h/a) of 0.17. Taking into account the orientation 
of the cluster in the tidal field of the Galaxy, Raboud & 
Mermilliod (1998) derive an intrinsic ellipticity of almost 0.3 
(see also van Lecuwcn ct al. 1986), comparable to what we 
find in our models (see Figure 12) by comparing the distance 
to the Jacobi surface along the z axis with the distance to 
Li: e ~ 1 - rz/rj. 



4.2 Praesepe 

Figure 17 shows the A/_R-magnitudc global luminosity func- 
tion for the W6 clusters at birth and at 800 Myr and com- 
pares that with the observed luminosity function for Prae- 
sepe reported by Hambly et al. (1995a, 1995b), which is 
shown as filled circles with error bars. 

The two 800 Myoar old luminosity functions are taken 
from the stars within the projected half mass radius and 
stars farther away. The outer luminosity function (dashes) 
fits better to the bright stars where the inner luminosity 
function fits better to the dimmer stars. Again we can argue 
that omitting the brightest stars from the sample provides 
a better fit to the observations, in which case the luminosity 



100 200 400 600 800 
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Figure 17. Cumulative luminosity function (Mjj-magnitude) of 
the Pra«sepe cluster. The filled circles with error bars give the 
observed luminosity function within the half-mass radius (Hambly 
at al. 1995a, 1995b), to which 170 stars with Mr < 4.5 were 

added. The model W6 luminosity function for all stars at zero 
age is shown as a dotted line. The solid and dahsed lines give 
the luminosity functions for models W6 at 800 Mycar within and 
outside the projected (on the y-z plane) half mass radii. 

function of the inner half of the cluster provides a better 
comparison. 

4.3 The Hyades 

Figure 18 compares the M/-band luminosity functions for 
the stars and binaries of several models, at birth (dotted 
lino) and at 600 Myr (other Hues), with the observed lu- 
minosity function the Hyades (Reid & Hawley 1999). Reid 
& Hawley observed the entire cluster and their luminosity 
function reportedly extends down to the hydrogen-burning 
limit. The W()=4 data for the entire cluster (solid) and the 
Wo =6 model for the stars within the half mass radius (dash- 
3dot) does not fit as well as the data for model Wo=4 for 
stars within the inner half mass radius (dashes). Our models 
W6 are somewhat farther from the Galactic center than is 
the Hyades, which would tend to suppress mass segregation 
somewhat by increasing the cluster tidal radius. The Wo=4 
models do excatly the opposite. Based on these arguments 
we conclude that Hyades is somewhat more mass segregated 
that our models predict and some degree of promiridal mass 
segregation seems to be required. 

Oort (1979) compared observations of Hyades with 
Aarseth's (1973, 1975) Af-body calculations and concluded 
that the outer 4 pc of the Hyades cluster are more strongly 
flattened, with e ~ 0.5, than the N-body models implied. 
The orientation of the Hyades in the Galcixy relative to the 
position of the sun then implies that the intrinsic flattening 
is even greater. Our calculations do not support Oort's con- 
clusion, and a flattening of e = 0.5 is quite consistent with 
our A^-body models. The reason for the discrepancy between 
our results and the conclusion of Oort is based on Aarseths' 




S 10 IS 

M, 

Figure 18. Cumulative M/-magnitude luminosity function (filled 
circles with error bars) of the Hyades star cluster (Reid & Hawley 
1999). The dotted line is the initial luminosity function for all 
stars in the model calculations. The solid line and dashed lines 

give the luminosity function for model W4 at an age of 600 Mcars 
for stars within the projected (on the y-z plane) tidal radius 
(solid) and within the half-mass radius (dahes). The dash-3dotted 
line gives the same data as the dashed line but then for model 
W6. 

models which were computed with a very small number of 
stars. The flattening of the cluster in the tidal fleld of the 
Galaxy, however, becomes more apparent towards the clus- 
ters' tidal radius which has smallest stellar density. Calcu- 
lations which are performed with a limited number of stars 
^ 500 will hardly show the falttening in the tidal field. 

4.4 NGC 3680 

Figure 19 compares the observed Mv-band luminosity func- 
tion of the old open cluster NGC 3680 with our model lumi- 
nosity functions. The observed luminosity function is poorly 
reproduced by our cluster models. However, if we remove 
the least luminous stars (those with V > 11.5), the 1.4 Gyr 
model flts the observed luminosity function fairly well. The 
imposed lower limit is rather arbitrary, it suggests that the 
observations may not properly correct for the faintest stars, 
or they may simply be absent from the data. Mass segrega- 
tion generally causes the lightest stars to escape from the 
cluster which, in time, leads to an over abundancy of mas- 
sive stars. The observed clusters, however, seem to have too 
few high mass stars and also too few low mass stars, which 
is hard to understand from a dynamically point of view. We 
therefore argue that in this case the lack of low mass stars 
is an observeational selection effect. 

4.5 Isophotes 

Figure 20 shows isophotes of the clusters NGC 2287, 
NGC 2516 and NGC 3680. NGC 3680 is most strongly flat- 
tened (e ~ 0.23), the other two are more circular in appear- 
ance; NGC 2287 has e - 0.05 and NGC 2516 has e ~ 0.14. 
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Figure 19. The filled circles with error bars show the observed 
cumulative luminosity function, corrected for field stars, of the 
star cluster NGC3680 (Hawley et al. 1999). The dotted line gives 
the initial luminosity function for all models, all other lines show 

luminosity functions at 1400 Mycar. The solid and dashed lines 
give the luminosity function for all stars within the half mass 
radius of model W4 and W6, respectively. The dash-3dotted line 
shows only the stars with My < 11.5 from the dashed line (within 
the half- mass radius of model W6 at 1400 Myear) . 



We measured these ellipticities for the inner 4pc for each 
cluster. 

The faintest stars in Figure 20 are 13th magnitude 
for NGC2287 and 15th magnitude for NGC2516 and 
NGC 3680. Taking the distances to these clusters into ac- 
count, it is clear that only the top end of the main sequences 
are included in the isophotes. Since mass segregation causes 
these high-mass stars to be more centrally concentrated than 
lower-mass stars, it is likely that we see only the inner, 
roughly spherical, regions of NGC 2287 and NGC 2516 (see 
Figure 12). NGC 3680 appears much somewhat rectabgular 
on the image and we are tempted to explain this shape as a 
result of the limited imaging; is the cluster is truncated by 
the edge of the field of view?. In this case it is not surprising 
that this cluster appears somewhat flattened. If one looks at 
the inner density contours the cluster looks much more cir- 
cular. Note also that the tidal radii for these clusters listed 
in Tab. 1, 13, 6.3 and 4.3 pc for NGC 2516, NGC 2287 and 
NGC 3680 respectively, seem inconsistent with the above 
pictures. The isophotes of NGC 2516 and NGC 2287 repre- 
sent only the central portion, where for NGC 3680 appears 
much bigger than the tidal radius we derived in Tab. 1. The 
table possibly underestimates the tidal radius of NGC 3680. 



5 COMPARISON WITH OTHER WORK 

Table 10 compares the evaporation times of our model cal- 
culations with previously reported results. We discuss each 
in turn. Table 10 lists the initial conditions of the models 
with which we compare our own. 

Terlevich (1987) performed direct iV-body calculations 
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Figure 20. Isophotes in My of the cluster NGC 2516 (top panel) 
NGC 2287 (bottom left), and NGC 3680 (bottom right). The 
adopted distance moduli are 9.1, 8.2, and 9.95, corresponding 
to distances of 655 pc, 373 pc and 735 pc for the top and bot- 
tom left and right panels, respectively. The (sub)giants are plot- 
ted as filled circles; the remaining stars are plotted as isophotes. 
The central luminosity densities are -2.0 mag pc^'^, -2.2magpc~-^ 
and -0.54 mag pc~^ for the top, bottom left and right panels, re- 
spectively. Contours are plotted at constant intervals of 1.25 mag. 
Stars are assigned a Gaussian point spread function with a disper- 
sion of 0.06 pc (0.1 pc for the bottom figure), the figures contain 
594, 178 and 738 stars for the top and bottom left and right 
panels, respectively. 
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Table 10. Overview of the model features, initial conditions and lifetimes of our model calculations and those reported by other workers. 
The Galactic tidal field may be modeled as that of a self-consistent disc (Disc) or a point mass (PM); initial density profiles are an 
anisotropic King model ( AKing) , an Isothermal sphere (Iso) , or a Plummer sphere (Plummer) the latter two with a cut off. Stellar and 
binary evolution are indicated by + (included) or - (for neglected) . Note: Terlevich and de la Puente Marcos work in terms of the initial 
virial radius, which is typically about 20% bigger than the half-mass radius. 





Galaxy 


Density 


Stellar Binary 


N 


(m) 


fhin 


''tide 


''hm 




^end 


^end/^rlx 




model 


profile 


evolution 




[Mq] 




[pc] 




[Myr] 




W4 


Disc 


AKing 


+ + 


3072 


0.56 


50 


6.3 


2.14 


109 


1200 


11.0 


W6 


Disc 


AKing 


+ + 


3072 


0.54 


50 


12.1 


2.00 


102 


~1600 


15.7 


Terlevich 


Disc 


Iso 


+ 


1000 


0.50 





12.1 


2.00 


71 


~1400 


19.7 


McMillan & Hut 


PM 


Plummer 




2048 


1 


10 


16.5 


2.06 


86 


2500 


29.0 


Kroupa 


PM 


Plummer 


+ 


400 


0.32 


100 


6.5 


2.53 


77 


693 


9.0 


de la Puente Marcos 


Disc 


Iso 


+ 


750 


0.60 


33 


9.9 


2.47 


67 


1061 


15.8 



with up to 1000 stars, including a power-law stellar mass 
function, mass loss from single-star evolution, and the tidal 
field of the Galaxy. The implementation of the Galactic tidal 
field and the evolution of single stars were somewhat sim- 
ilar to those presented in this paper. Her model XII had 
initial conditions most similar to our own, although some 
significant differences exist. Terlevich's models started with 
spherical distribution of stars with density proportional to 
virial radius 2 pc (half-mass radius of about 1.6pc), 
and located at a distance of 10 kpc from the Galactic center. 
The resulting half-mass crossing time was 5.2 Myr, compa- 
rable to the ~ 4.1 Myr in our models. Terlevich's model XII 
initially consisted of 1000 stars drawn from a Salpeter mass 
function, with a mean mass of 0.5 Mq. The initial virial ra- 
dius was Q = -Ekin/|Spot| = 0.25, less than equilibrium value 
of 0.5. Even though the model started with fewer stars, the 
cool initial conditions mean that the initial relaxation time 
was comparable to that in our own models. The half-life 
of model XII was 770 Myr. The run was terminated at 1 
Gyr, by which time about 70% of the mass had been lost. 
We estimate that the cluster would have dissolved in about 
1.4 Gyr. 

The half-mass lifetimes of our W6 models are all around 
800 Myr, similar to that of Terlevich's model XII. This is 
somewhat surprising, as we might have expected that our 
W6 models, with more stars, would live longer than the com- 
parable models of Terlevich. The reason for this discrepancy 
cannot be attributed to the large fraction of primordial bi- 
naries in our calculations, as binaries do not dramatically 
affect the evaporation rate of the cluster (Figure 1). The 
difference may possibly bo due to the lower high-mass cutoff 
in Terlevich's initial models. Although her model XII was 
initially less concentrated than our W6 models, it did reach 
core collapse. In general, more concentrated models tend to 
live considerably longer than shallower models (Takahashi 
& Portegies Zwart 2000). 

The evolution to core collapse in Terlevich's model XII 
was aided by the absence of stars with masses greater than 
10 M0. By the time the turnoff mass has dropped below 
10 Mq (after about 22 Myr) our W6 models have lost be- 



tween 4% and 9% of their mass due to stellar evolution 
alone. The loss of even such a small mass fraction may 
have dramatic consequences for the further evolution of the 
stellar system, as this mass is lost from the most massive 
stars, which reside deep inside the cluster potential well. The 
shorter lifetime of clusters with large populations of massive 
stars is demonstrated by Terlevich's model XV, where the 
initial mean mass was 7.4 MQand the cluster did indeed dis- 
solve much more rapidly. 

The Galactic tidal field produced a similar flattening 
effect in Terlevich's clusters as in our own (compare her 
Figure 7 with our Figure 12). 

The simulations reported by McMillan and Hut (1994) 
included up to 2048 equal mass stars, including up to 20% 
rather soft primordial binaries, and incorporated the Galac- 
tic tidal field, modeled as the field of a distant point mass. 
However, they excluded stellar evolution and iieiicc any pos- 
sibility of stellar mass loss. In the absence of a physical 
time scale associated with stellar evolution, they presented 
their results in units of the initial half-mass relaxation time. 
Our W6 models have half-lives of about 6 initial relaxation 
times, much shorter than the ~ 29 initial relaxation times 
for the most comparable McMillan & Hut models. Also, as 
discussed previously, all the McMillan & Hut models expe- 
rienced core collapse, which is absent in our simulations. 

The main reasons for these differences are the effects of 
stellar mass loss and the presence of a stellar mass function 
in the present studies. In our models, core collapse is arrested 
by mass loss. In addition, the McMillan & Hut models all 
started off well inside their tidal radii, significantly increas- 
ing their lifetimes. 

In a continuing effort to understand the evolution of 
young open star clusters, Kroupa (1995a; 1995b; 1995c) per- 
formed Af-body calculations with up to 400 stars, all of them 
members of primordial binaries. He adopted a distribution in 
orbital separation flat in log a, but selected a between ~ 360 
and ~ 3.6 10^ R©; the binaries in his calculations were thus 
on average much wider than in our calculations. His models 
included the Galactic tidal field and stellar mass loss, but 
neglected binary evolution. 
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Kroupa's model with the highest initial relaxation time 
started from a Plummer sphere with a half-mass radius of 
2.5 pc and a crossing time of ~ 20 Myr. It dissolved in about 
700 Myr. Scaled to the initial relaxation time, this is some- 
what faster than our W4 models. The main reasons for this 
more rapid evaporation are most probably the shallower ini- 
tial density profile in his models and the cluster's smaller 
distance from the Galactic center. 

We compare the binary and triple properties of our 
models with those of McMillan & Hut and Kroupa in more 
detail in Paper IVb. 

De la Puente Marcos (1997) studied the effect of the ini- 
tial mass function on the dynamical evolution of open star 
clusters, including both stellar mass loss and the tidal field of 
the Galaxy. His calculations were limited to 750 stars and in- 
cluded 33% rather wide primordial binaries, all with a mass 
ratio of 0.5, in which (binary) evolutionary effects were ne- 
glected. His model XX used Scale's (1986) mass function and 
had an initial virial radius of 2.47 pc. The initial relaxation 
time for this model was about 67 Myr; the crossing time was 
7.5 Myr. This model dissolved in about 1 Gyr, slightly faster 
than our models. However, scaled to the initial relaxation 
time, this result is consistent with the dissolution time of 
our W6 models (see Table 10). 



6 SUMMARY AND DISCUSSION 

The aim of our simulations is to study the evolution of open 
star clusters such as the Pleiades, Praesepe and the Hyades. 
These clusters differ significantly in age, but have compara- 
ble physical characteristics, stellar membership, total mag- 
nitude and internal velocity dispersions. Our TV-body calcu- 
lations incorporate, in a fully self-consistent fashion, mass 
loss from single stars, binary evolution, dynamical encoun- 
ters among single stars and binaries and the effect of the 
Galactic tidal field. 

We have compared the luminosity functions, isochrones 
and projected luminosity profiles of our models with obser- 
vations. For some clusters, it is hard to find a good match 
between model and observed luminosity functions. Mass seg- 
regation and observational limitations significantly reduce 
our ability to find a match, and restrict our understanding 
of the differences we see. However, we find that for models for 
which we compared the luminosity functions with Preiesepe 
and the Hyades, the observations show evidence for signif- 
icantly more mass segregation than is seen in the models. 
We conclude that these clusters may have been born some- 
what mass segregated. Alternatively these clusters may have 
started out somewhat more massive than assumed here, but 
with a shallower density profile. The selective evaporation 
of lower mass stars then results in a "dynamically old" ap- 
pearance (see also Takahashi & Portegies Zwart 2000). 

The global luminosity function of Praesepe and its 
degree of mass segregation suggest an ago greater than 
800 Myr, which is at the high end of the observed range. 
The luminosity function of NGC 2287 is consistent with the 



observed age of 150-200 Myr. Our age estimates for the 
Pleiades and Hyades, based on the structure and dynamical 
state of these clustusters, are consistent with ages derived 
from isochrone fitting. 



6.1 Mass segregation 

The first effects of mass segregation in our models are dis- 
cernible in the cluster core after only a few million years, a 
small fraction of an initial half-mass relaxation time. After 
about a relaxation time, mass segregation becomes measur- 
able in the cluster outskirts. The luminosity function of the 
inner parts of the cluster provides a useful tool for studying 
mass segregation, although one has to select regions which 
arc well separated in radius in order to make the effect visi- 
ble. As expected, giants and binary stars are most strongly 
affected by mass segregation and it is easiest to identify the 
effect by comparing the radial distribution of giants with 
that of the lower mass main-sequence stars. 

One important effect of mass segregation is that older 
clusters become rich in white dwarfs and giants relative to 
the Galactic field. These stars may be single or members of 
binary systems. The main reason for the overabundance of 
giants and white dwarfs in clusters is the the depletion of low 
mciss main-sequence stars by evaporation. This fiattening of 
the mass function due to mass segregation has also been 
studied by Takahashi & Portegies Zwart (2000) for more 
massive clusters. They also find that clusters which arc close 
to disruption are rich in compact objects and giant stars. 

The W4 models have shallower initial potentials, lie 
closer to the Galactic center, and are more strongly affected 
by mass segregation because of their more rapid evapora- 
tion. We suggest that the best place to look for evidence 
of mass segregation is in star clusters with larger half-light 
radii, which have shallower potentials. Since mass segrega- 
tion manifests itself more clearly in dynamically evolved sys- 
tems, it is also better to look at older clusters. A relatively 
old cluster (age ^ 500 Myr) with a relatively small mass 
(Mtot ~ 500 M©) would be ideal. 



6.2 Age estimates 

The "dynamical ages"^ of our model clusters are often in- 
consistent with ages determined by isochrone fitting. The in- 
stantaneous relaxation time is a poor estimator of a cluster's 
dynamical age. The cores of star clusters lose their memory 
of the initial relaxation time within a few million years, well 
within an initial half mass relaxation time, but the time scale 
for global cluster amnesia to set in is far larger than the ini- 
tial relaxation time. The half-mass relaxation time tends to 
increase by a factor of two or three, reaching a maximum 



^ We define the dynamical age of a cluster as its lifetime ex- 
pressed in units of the initial relaxation time. 
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near the cluster's half-life and decreasing thereafter. The in- 
crease is caused by the internal heating of the cluster; the 
decrease mainly by loss of stars. 



6.3 Escaping stars 

Stars tend to escape in the direction of the Li and L2 La- 
grange points of the Galaxy-cluster system. The velocities 
of the stars at these points arc highly anisotropic and, as ex- 
pected, pointing mostly radially away from the cluster cen- 
ter. The velocities of the escaping stars arc comparable to 
the cluster velocity dispersion; very few stars are ejected 
with high velocities following a strong dynamical encounter. 

Neutron stars arc ejected isotropically and with much 
higher velocities, owing to the asymmetric kicks they receive 
during the supernovae that create them. Binaries contain- 
ing the progenitors of neutron stars are generally disrupted 
by the first supernova explosion; in all the calculations pre- 
sented in this paper, only one binary survived the first su- 
pernova. However, the existence of that single binary does 
suggest that it may be possible to form X-ray binaries in 
open clusters. Such binaries are only expected to exist in 
star clusters which are younger than ~ 45 Myr (the turnoff 
age of a 7Mq star), because mass loss and the velocity kick 
imparted to the binary causes it to escape (see paper IVb). 

Black holes are more easily retained by the cluster, but 
are very rare due to their high progenitor mass and the 
steepness of the Scalo initial mass function. 



6.4 Tidal flattening 

The tidal field of the Galaxy fiattens the cluster significantly 
in the z direction, and to a lessor extent along the y axis. A 
cluster which is spherical at birth develops this flattening in 
its outer regions within a few crossing times; the inner parts 
remain fairly spherical. All our models show this flattening, 
but its observability from Earth depends on the orientation 
of the cluster in the Galactic plane. 

EUipticities reported for the Pleiades (e ~ 0.3; van 
Leeuwen et al. 1986) and Hyades (e ^ 0.5 in the outer 
regions; Oort 1979) arc consistent with our model calcu- 
lations. The available data for NGC 2287 and NGC 2516 do 
not show signiflcant flattening. However, these data contain 
only stars from the innermost 3 pc, which are affected least 
by the Galactic tidal field. NGC 3680 appears more square 
than elliptical. This may be caused by the small field of view 
of the telescope, which could cause a star cluster to take on 
the shape of the CCD frame (see Fig. 20) . 

The flattening persists during mass segregation, in the 
sense that the mass-segregated Lagrangian radii are ellip- 
soids. Projection of the cluster onto the background may 
therefore decrease the observed mass segregation. 



6.5 Core collapse 

Our models experience rather shallow core collapse during 
their early evolution. The clusters then expand more or less 
homologously, preserving their initial density profiles. The 
expansion is driven by stellar mass loss and, to a lesser ex- 
tent, by binary heating; dynamical models without stellar 
mass loss but which include a tidal field and primordial bi- 
naries do show core collapse (McMillan & Hut 1994). Once 
the cluster has lost a considerable fraction of its stars the 
system shrinks again. The remnant with a few remaining 
stars may become quite compact before it dissolves, but we 
find no evidence for late core collapse before complete dis- 
ruption. 



6.6 Giants and white dweirfs in open clusters 

A cluster's single-star population is not noticeably affected 

by cluster dynamics until the age of the system exceeds 
^ 2 initial half-mass relaxation times. However, the binary 
populations are measurably influenced by dynamics oven at 
early times (see Paper IVb) . At later times {t > 400 Myr) , 
our model clusters tend to become rich in giants and white 
dwarfs. Small-number statistics on the giants limit the de- 
gree to which we can quantify this statement , but the white- 
dwarf populations in our models increased by factors of 1.3 
to 4.6 (for model W6 and W4, resepctively) relative to what 
one would expect for a non-dynamically evolving population 
of single stars and binaries. Most of the excess white dwarfs 
arc members of binary systems. 

Table 1 presents an overview of the numbers of white 
dwarfs observed in the clusters studied in this paper. Ex- 
plaining the number of white dwarfs in the Hyades is a long- 
standing problem, starting with discussions in the mid-1970s 
by Tinsley (1974) and van den Heuvel 1975), and continu- 
ing into the 1990s (Eggen 1993; Weidemann 1993). All these 
papers conclude that the observed number of white dwaxfs 
is too small, by about a factor of three, for the inferred clus- 
ter mass and age. The three loading explanations wore: (1) 
the upper mass limit for the production of white dwarfs 
may be as low as ~ 4M0 (Tinsley 1974), white dwarfs 
are born with a velocity kick as neutron stars do (Weide- 
mann et al. 1992) and (3) mass segregation selectively ejects 
white dwarfs (van den Heuvel 1975). By studying the white 
dwarfs in NGC 2516, Koester & Reimers (1996) firmly con- 
clude stars up to 8Mq can form white dwarfs, removing 
the flrst solution to this conundrum. However, our calcula- 
tions arc inconsistent with the idea that white dwarfs are 
preferentially ejected from star clusters. On the contrary, 
the white dwarfs in our simulations are more easily retained 
than main-sequence stars of the same mass, causing the older 
clusters to become relatively white-dwarf rich for their mass. 

We can study this problem further by comparing the 
ratio of the number of white dwarfs to the number of giants: 
f wd = A^wd/A^gs. For the first five open clusters in Tab. 1 this 

fraction ranges from /wd = 1 for NGC 2516 and the Pleiades 
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to /wd = 2.5 for the Hyades. For an evolving population of 

gs 

single stars without dynamics (see Tab. 2), we find that the 

ratio ranges from 1.1 at 100 Myr (the age of NGC 2516) to 
4.2 at 600 Myr (comparable to the age of Hyades). 

Binary evolution complicates the comparison due to an 
obvious selection effect — the giants can probably all be seen, 
but white dwarfs are easily hidden near a main-sequence or 
giant companion. We obtain estimates for /w<i in a field 

population with 50% primordial binaries by combining the 
single stars from Tab. 2 with the binaries from Tab. 8. We 
calculate an upper limit to /wd by accounting for all white 

gs 

dwarfs [counting (wd, wd) binaries as two objects and in- 
cluding (ms, wd) and (gs, wd) binaries]. A lower limit is 
obtained by counting (wd, wd) binaries as one object and 
excluding white-dwarf binaries containing a main-sequence 
or giant companion. For the field (no dynamics) population, 
we find /wd = 0.80-1.3 at 100 Myr, and /wd = 2.7-4.1 at 

600 Myr. Combining models W4 (Tab. 7) and W6 (Tab. 6) 
we obtain /wd = 1.2-1.7 at 100 Myr and /wd = 3.0-4.0 at 

gs gs 

600 Myr. 

The observed value of /wd ~ 1 at around 100 Myr seems 

somewhat low, but probably not inconsistent with our mod- 
els. However, the value of /wd ~ 2.5 of the Hyades cluster 

gs 

(at about 600 Myr) is smaller than our models predict, sug- 
gesting that a considerable fraction of the white dwarfs are 
hidden in binaries. 

The value of the fraction /wd does not seem to pose a 

gs 

serious problem to understand any of the clusters discussed 
here. But according to the evolution of the field stars and 
binaries (combining Tabs. 2 and 8), at an age of 600 Myr 
we expect a total of about 15 giants and 59 white dwarfs. 
The comparable models W4 and W6 contain, respectively, 
13 and 12 giants, and 48 and 57 white dwarfs at that age. 
Model W6 has lost about 40% of its mass and model W4 
has lost about 60%, yet the numbers of giants and white 
dwarfs have decreased by only 4% - 20%. Averaging over 
time, we find a decrease in the rmmber of giants and white 
dwarfs of about 2% per 100 Myr. Apparently, the dynamical 
evolution of these clusters has little effect on the number 
of giants and white dwarfs. (In fact, this was assumed by 
von Hippel [1998] in estimating the mass fractions of white 
dwarfs in open clusters.) Number counts of giants and/or 
white dwarfs may therefore provide a reasonable estimate of 
a cluster's initial mass. 

For a fifty-fifty mixture of single stars and binaries, 
meaning that 2/3 of the stars are binary components, the 
mean mass of a star is (m) = 0.46 at t = 100 Myr, 0.41 
at 600 Myr and almost constant ((m) ~ 0.40) thereafter. 
We use the numbers of giants to estimate the initial masses 
of the star clusters in Tab. 1, because the giants are least 
plagued by selection eflfects (although their small numbers 
significantly limit the accuracy of our estimates). The num- 
ber of giants per Ik stars is 2.7 at 100 Myr (see Tabs. 2 and 
8), and rises rapidly to about 6.5 at 400 Myr, after which 



the specific number of giants remains roughly constant. For 
an open cluster older than ~ 400 Myr, we thus estimate its 
intial mass via 

Mo=65MeiV,.(l + ^^). (1) 

For younger clusters, the factor 65 (= 1024 x 0.41/6.5) is 

larger 170 at 100 Myr and ~ 100 at 200 Myr. 

We apply this method to the clusters from Tab. 1 with 
more than 5 giants and obtain the following birth masses: 
830 Mo for NGC 2287, 3100 Mq for NGC 2660, and 1400 Mq 
for NGC 3680. These mass estimates seem reasonable. For 
NGC 2534, Praesepe and the Hyades, the mass estimates 
are 690 M©, 370 M0 and 290 M©, respectively, consider- 
ably smaller than the observed masses of these clusters (see 
Tab. 1). 

The initial mass estimate for these clusters increases 
proportional to the number of giants. The ratio /wd does 

gs 

not pose a serious problem and therefore, instead of too 
few white dwarfs, Hyades may have too few giants. Where 
white dwarfs can be hidden easily in binaries, giants are 
not that easy to hide. One way to decrease the number of 
giants is by binary activity. The number of giants can be 
decreased when subgiants are stripped in a phase of mass 
transfer before they reach the horizontal branch, where they 
spend most of their time. In order to reduce the rmmber of 
giants in the fashion we require most binaries to be born 
with short orbital periods ( ^ 100 days). Alternatively the 
giant lifetimes adopted in our models may be too long. 

6.7 Notes on individual clusters 

NGC 2516 deserves much more study, as its dynamical pa- 
rameters (total mass, half-mass and core radii) are poorly 
known. 

The Pleiades cluster has been quite thoroughly stud- 
ied in searches for brown dwarfs and planets. Its mass func- 
tion fits nice with the luminosity function from our models, 
with the exception that our models containt 00 many bright 

stars. 

NGC 2287 is quite poorly studied, and with the quoted 
values for the tidal radius it is has an extremely low mass for 
its age (see Table 1). The presence of 8 giants suggests that 
its mass must have been comparable to that of NGC 2516. 

Praesepe has been studied recently in great detail, and 
appears to fit well with our dynamical models. However, the 
cluster is rather shallow and may have been born somewhat 
more massive and less concentrated (Wo ^ 4) than our mod- 
els. This conclusion is based on the observed shallow den- 
sity profile and the high degree of mass segregation. There 
is some excess of stars with V 9-11 is unexplained by our 
models. The observed cluster seems to be very deficient in 
gaints. Based on the number of white dwarfs and the dynam- 
ical state of the cluster we sould expect at least 20 giants in 
this cluster, where only 5 are known. 

The Hyades cluster fits well with our models, indicat- 
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ing that it is possible to estimate the initial conditions for an 
observed star cluster rather accurately. The cluster does not 
appear to be deficient in white dwarfs if we compare them to 
the number of giants. However, if the mass of Hyades quoted 
in Table 1 is correct, the observed number of white dwarfs 
and giants seems too small by a factor of about three. 

NGC 3680 Fits well with our W6 models expect for 
the luminosity function, which is deficient of low mass stars. 
A possible solution may be that these stars were initially 
absent in the cluster or that the observations do not go faint 
enough to reveal the low mass stars. 

6.8 Comparison with other work 

Our models evaporate on time scales generally consistent 
with dissolution times reported in previous calculations. The 
models of Terlevich (1987) and de la Fuente Marcos (1997) 
compare well with the evaporation times of our W6 models. 
Terlevich's models evolve somewhat more slowly due to the 
lack of massive stars and the rather cool initial conditions, 
which drive the cluster to core collapse and therefore ex- 
tend its lifetime somewhat. The models of McMillan & Hut 
(1994) dissolve more slowly than ours. This discrepancy can 
be completely explained by the absence of stellar evolution 
mass loss and binary evolution in their models, along with 
the small size of those models relative to the clusters' tidal 
radii in the Galactic potential. 

The only real discrepancy is in the work of Kroupa 
(1995a), whose models dissolve somewhat more rapidly than 
ours. Possibly the small numbers of stars and the large 
binary fractions drive a more rapid evaporation than one 
might naively expect. The evaporation rate of star clusters 
is known to depend on total the number of stars (Heggie et 
al. 1997; Portegies Zwart et al. 1998). In Section 1 we argue 
that the presence of primordial binaries has little effect on 
the cluster lifetime. It is, however, not clear how this trend 
propagates in the scaling of cluster lifetimes with respect to 
the number of stars. 



APPENDICES 



APPENDIX A: TERMINOLOGY 

Throughout the paper (and in future papers in this series) 
we will use consistent nomenclature. Some of these terms 
arc rather confusing and have been used by different authors 
in the past with slightly different meanings. For clarity we 
present here a short glossary of terms. 

Binary fraction: thoughout this paper we define the 
binary fraction as: 



fhin — 



Afbin 



A^'sing + A^bin ■ 



(Al) 



Here A'^sing and ATbin are the number of single stars and bi- 
naries. 

Cluster center: The density center of the cluster, as 

defined below. Alternative definitions may use the number 
density or the luminosity density, or the point where the 
density is greatest, and may include projection effects. 

Collision: A collision occurs when the distance be- 
tween two stars i and j becomes smaller than the sum of 
their effective radii: d < dcou{ri + rj), with dcoii = 1- The 
effective radii arc determined from detailed fluid-dynamical 
calculations. 

Core radius: The weighted average distance of all stars 
from the density center. Casertano & Hut (1985): originally 
used a weighting proportional to the local density. However, 
in practice this definition is unsuitable for clusters with near- 
isothermal density profiles (p ^ Following Aarseth 
(1986), we adopt the modified definition 



(A2) 



where pj is given by Eq. A6. Note that this definition of 
the core does not necessarily have any simple relation to 
the "core radius" normally quoted by observers, nor to the 

"dynamical" core radius Vc = \J^^^f-, where pc and {v'^)c 
are, respectively, the cluster's central density and velocity 
dispersion (Binney & Tremaine 1987). 

Crossing time: The time taken by a star with velocity 
equal to the velocity dispersion v to cross the virial radius 
r of the stellar system tc = r/v which for practical reasons 
is written as 

^'-[gmJ ■ 

In more convenient units, we write the half-mass crossing 
time as 



(A3) 



(A4) 
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Density center The density weighted average of the 
positions of all stars (von Hoerner 1960, 1963): 



VH 
^dens 



(i) 



j 



(A5) 



In these expressions, p^*' is the density estimator of order 
j around the i-th particle, with position vector rj. For any 
star i we define the local density within the volume Vj of 
the sphere containing the j nearest neighbors {11,12, ■■■ 

of i as 



Vi 



where Vi 



3 3 l'» 



(A6) 



and the sum over masses excludes 



the masses of both stars i and ij (see Casertano & Hut 1985). 
We take j = 12. 

Escaper: A star which is not bound to the cluster, 
i.e. whose energy exceeds the energy at the cluster's Jacobi 
surface. An escaper may lie within the Jacobi surface. 

Half mass radius: The radius of the sphere, centered 
on the cluster density center, that contains half of the total 
cluster mass (as defined in the text). It is not always clear 
which stars to include in determining this radius, as the 
cluster is generally flattened in the Galactic tidal field. 

Hcird binary: A binary whose binding energy exceeds 
the mean stellar energy in the cluster (Heggie 1975). A bi- 
nary is hard if its semi major axis exceeds 



GMm{M + m + ma) 



(A7) 



(M + m)m3V^ 

Jacobi radius: The distance from the cluster center 

to the Li and L2 Lagrange points- the maximum distance 
from the cluster center to the cluster Jacobi surface. 

Jacobi surface: The cluster's "Roche lobe" in the tidal 
field of the Galaxy. Consider a star moving with Jacobi inte- 
gral Ej = ^v'^ + 4>ef f{v) in the rotating frame of reference in 
which our simulations are performed, where v is velocity and 
the effective potential 4>ef f includes the cluster's self-gravity, 
the tidal field of the Galaxy, and the centrifugal force in the 
rotating frame. The zero-velocity surface for that value of Ej 
is defined by u = 0, so 0e//(r) — Ej (Binney & Tremaine 
1987). The Jacobi surface for the cluster is defined to be the 
last closed zero-velocity surface that contains the cluster — 
that is, the surface passing through the cluster's Li and L2 
Lagrange points. With the conventions adopted in kira, the 
Lagrange points are located along the x-axis, the cluster or- 
bits in the x-y plane, and the Jacobi surface is elongated 
in the x-direction. The three coordinate axes intersect the 
Jacobi surface at distances Vx = rj (the Jacobi radius), r^, 
and rz from the cluster center. 

Member: A star (or multiple system) which is bound 
to the cluster. The energy of such a star is less than the en- 
ergy at the cluster Jacobi surface. A member may lie outside 
the Jacobi surface. 

Primary: The more massive of the two stars in a binary 
system. Denoted with M. 



Relaxation time: We use Spitzer's (1987) definition 
of the half mass relaxation time: 



GM 



1/2 



N 



5 log A ■ 



(A8) 



Here A ~ 0.4A'^ is the Coulomb logarithm. (In the presence 
of a realistic mass function, A ^ O.IA'^ may be more appropri- 
ate; Farouki & Salpctcr 1982; 1994; Smith, 1992; Fukushige 
& Heggie 1999). In convenient units this may be written 



trlx = 2.05 



[Me 



1/2 



[pc] 



3/2 



N 
log A 



[Myr]. 



(A9) 



Note that "iV" here is the number of bound objects in the 
star cluster, and is smaller than the total number of stars if 
the cluster contains binaries. 

Secondary: The less massive of the two stars in a bi- 
nary system. Denoted with m. 

Tidal radius: Since clusters are somewhat elongated 
in the Galactic tidal field, the tidal "radius" is not well de- 
fined. For definiteness, we take the tidal radius to be the 
Jacobi radius of the cluster. For a disk field described by 
Oort (1927) constants A and B, we have 

GMtot 



[pc] 



(AlO) 



4:A{A - B) 

In the solar neighborhood, A = 15 km s~^kpc~^, and B = 
-12 km s"^kpc"^ 

Unperturbed binary: A binary for which the dimen- 
sionlcss perturbation due to its neighbors is less than some 
critical value (~ lO"'', typically). Dynamically, unperturbed 
binaries are treated in the point-mass approximation, as seen 
by the rest of the system. Only unperturbed binaries can be 
treated using the SeBa binary evolution module described 
in Appendix B. Perturbed binaries are treated as two single 
stars; mass transfer and tidal circularization are currently 
not handled in perturbed binaries. 

Virial radius: A characteristic length scale for the sys- 
tem, defined by 



-2U 



(All) 



where Afr/ is the total potential energy of the system (includ- 
ing the tidal potential). For an isolated equal-mass system, 
Tvir is the harmonic mean of the particle separations (Henon 
1972). 



© 0000 RAS, MNRAS 000, 000-000 



24 Portegies Zwart et al. 



APPENDIX B: STARLAB 

The simulations described in this series of papers are carried 
out within the "Starlab" software environment, version 3.5. 
Starlab is a software package for simulating the evolution 
of dense stellar systems and analyzing the resultant data. It 
consists of a collection of loosely coupled programs ( "tools" ) 
linked at the level of the UNIX operating system. The tools 
share a common data structure and can be combined in arbi- 
trarily complex ways to study the dynamics of star clusters 
and galactic nuclei. The main components of Starlab used 
in this work are kira, the A'^body integrator, and SeBa, a 
stellar and binary evolution package. The Starlab system is 
described in detail in http : / /www . sns . ias . eclu/~starlab. 

Bl Kira 

The iV-body integrator kira is the laxgest single program 
within Starlab. Its basic function is to take an input W-body 
system and evolve it forward for a spccifiod period of time, 
producing snapshot and other diagnostic output at regular 
intervals. In addition to strictly dynamical evolution of stars 
and multiple stellar systems, kira also incorporates stellar 
and binary evolution (via the SeBa subpackage), and the 
possible influence of an external ("tidal") gravitational field. 
The program is designed to take advantage of the "GRAPE- 
4" special-purpose processor (Makino et al. 1997), if avail- 
able, although GRAPE is not required for its operation. 

Bl.l The Integrator 

Particle motion is followed using a fourth-order, block- 
timestop (McMillan 1986) "Hermitc" predictor-corrector 
scheme (Makino and Aarseth 1992). Briefly, during a time 
step 5t, particle positions x and velocities v are first pre- 
dicted using the known acceleration a and "jerk" j (the time 
derivative of the acceleration): 

Xj, = yL + v6t+ ^aSt^ + ij(5f^ , (Bl) 
Vp = w + a5t+^j5t\ (B2) 

The acceleration ap and jerk jp are then computed at the 

predicted time using Xj, and Vp, and the motion is cor- 
rected using the additional derivative information thereby 



obtained, 

k = ^a"5t^ = 2(a-ap)-F<5t(j-jp), (B3) 

I = ia"'5t' = -3(a-ap)-5t(2j-Kjp), (B4) 
to obtain the corrected position and velocity: 

X, = xp-h(^l+ik)<5^^ (B5) 

vc = vp + (il+ik)<5t. (B6) 



A single integration step in thus proceeds as follows: 

(i) Determine which stars are to be updated next. Each 
star has associated with it an individual time t, representing 
the time to which it was last advanced, and an individual 



timestep St. The list of stars to be integrated consists of 
those with the least value of i -|- St. Time steps are con- 
strained to be powers of 2, allowing "blocks" of many stars 
to be advanced simultaneously. 

(ii) Before the step is actually taken, check for 

(a) termination of the run 

(b) escaper removal 

(c) system reinitialization 

(d) diagnostic ("log") output, which includes 

- information on bulk parameters of the system: 
total mass, energy, momentum, anisotropy, etc. 

- technical information on CPU time, timestep 
distribution, etc. 

- detailed information on the cluster mass distri- 
bution: core properties, Lagrangian radii, etc. 

- stellar mass distribution and anisotropy by Lar 
grangian zone 

- luminosity profile, and mass and luminosity 
functions 

- cluster stellar content (by spectral type and lu- 
minosity class) 

- detailed dynamical and physical data on all bi- 
nary systems. 

(e) snapshot output, for restart and display 

(iii) Perform low-order prediction of all particles to the 
new time. This operation may be performed on the GRAPE, 
if present. 

(iv) Recompute the acceleration and jerk on all stars in 
the current block (using the GRAPE, if available), and cor- 
rect their positions and velocities for fourth-order accuracy. 

(v) Check for and initiate unperturbed motion. 

(vi) Check for collisions and mergers. 

(vii) Check for tree reorganization (see below). 

(viii) Check for and apply stellar and/or binary evolution 
(§B.2), and correct the dynamics as necessary. 

B1.2 Tree Structure 

An iV-body system in Staxlab is represented as a linked- 
list structure, in the form of a mainly "flat" tree having 
individual stars as leaves. The tree is flat in the sense that 
single stars (i.e. stars that are not members of any multiple 
system) are all represented as top-level nodes, having the 
root node (the system center of mass) as parent. Binary, 
triple, and more complex multiple systems are represented 
as binary trees below their top-level center of mass nodes. 
The tree structure determines both how node dynamics is 
implemented and how the long-range gravitational force is 
computed. 

Each parent node contains "local" information about its 

dynamics — mass, position, velocity, etc. relative to its par- 
ent node. The leaves contain additional information about 
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stellar properties — effective radius, luminosity, temperature, 
etc. The parent node of a unperturbed binary also contains 
information on the binary parameters — semi-major axis, ec- 
centricity, mean anomaly, etc. The motion of every node 
relative to its parent node is followed using the Hcrmitc 
predictor-corrector scheme just described. The use of rel- 
ative coordinates at every level ensures that high numerical 
precision is maintained at all times, even during very close 
encounters. 

The tree evolves dynamically according to simple 
heuristic rules: particles that approach "too close" to one 
another are combined into a center of mass and binary node; 
and when a node becomes "too large" it is split into its bi- 
nary components. These rules apply at all levels of the tree 
structure, allowing arbitrarily complex systems to be fol- 
lowed. In practice, the term "too close" is taken to mean that 
two stars (1 and 2) approach within the "close-encounter dis- 
tance" -Rciosc ~ rvii{mi +m2)/2Mtot, the impact parameter 
that would lead to a 90° deflection if both bodies moved 
at typical stellar speeds. "Too large" means that a node's 
diameter exceeds 2.5i?ciose- 



B1.3 Binaries 

How the acceleration (and jerk) on a particle or node is com- 
puted depends on its location in the tree. Top-level nodes 
feel the force due to all other top-level nodes in the system. 
Forces are computed using direct summation over all other 
particles in the system; no tree or neighbor-list constructs 
are used. (This procedure is designed specifically to allow 
efficient computation of these forces using GRAPE hard- 
ware, if available.) Nearby binary and multiple systems are 
resolved into their components, as necessary. 

The internal motion of a binary component is naturally 
decomposed into two parts: (1) the dominant contribution 
due to its companion, and (2) the perturbative influence of 
the rest of the system. This decomposition is applied re- 
cursively, at all levels in a multiple system. Since the per- 
turbation drops off rapidly with distance from the binary 
center of mass, usually only a few near neighbors arc sig- 
nificant pcrturbcrs of even a moderately hard binary. These 
neighbors arc most efficiently handled by maintaining lists 
of perturbers for each binary. Perturber lists are recomputed 
at time the center of mass is integrated. 

A further cflicicncy measure is the imposition of un- 
perturbed motion for binaries whose perturbation falls be- 
low some specifled value for all or part of an orbit. Unper- 
turbed binaries may be followed analytically for many orbits 
as strictly two-body motion; they are also treated as point 
masses, from the point of view of their influence on other 
stars. The use of the unperturbed approximation near the 
periastron of eccentric orbits was a key element in our deci- 
sion not to use cumbersome regularization schemes for the 
computation of binary motion. 

Because unperturbed binaries are followed in steps that 
are integer multiples of the orbit period, we can relax the 



perturbation threshold for unperturbed motion relative to 
that for a perturbed step (since most of the perturbative 
effects of nearby stars are periodic). Perturbed binaries are 
resolved into their components, both for purposes of deter- 
mining their center of mass motion and for determining their 
effect on other stars. Unperturbed treatments of multiple 
systems are also used, based on empirical studies of the star 
bility of their internal motion. A hierarchical system is re- 
garded as stable if (a) the external perturbation is less than 
some threshold value, and (b) each component is stable (or 
single), by the same criterion. 

"Lightly perturbed" binaries, having external pertur- 
bations within a factor of ~10 of the unperturbed thresh- 
old, are treated using a variant of the method described 
by Mikkola & Aarseth (1998), in which the internal mo- 
tion of the binary is artificially slowed and the perturba- 
tion is increased by the same factor. Briefly, the result is 
that long-term secular trends in the binary orbital elements 
are properly reproduced, while periodic perturbative terms 
are amplifled; the latter effect is suppressed by following the 
"slow" motion over an integral number of orbits. Our "slow" 
binary treatment differs from that of Aarseth mainly in that 
it is not coupled to a regularization scheme it is applied di- 
rectly to the unregularized equations of motion. In addition, 
we apply pairwise corrections to forces between perturbers 
and the binary center of mass in order to avoid spurious high 
derivatives caused by the mismatch between the (slowed) in- 
ternal motion and the (normal) external interaction. 

El. 4 Tidal Field 

The standard form of the external (tidal) potential is 

4>e^t = ^{aix^ + asz^) . (B7) 

This expression includes contributions from both the Galac- 
tic tidal field and the centrifugal force in the cluster's rotat- 
ing frame of reference. The Galactic center is assumed to lie 
along the negative x axis and the rotation vector ft is in the 
z direction. (We assume motion in a circular orbit in the x-j/ 
plane around the Galactic center.) The equations of motion 
also include a Coriolis acceleration a.c — — 2f2 x v. Tidal 
and coriolis effects are applied to top-level nodes only — that 
is, we neglect the tidal effect of the Galaxy on a binary's 
internal motion. 

The values of ai, as, and fl depend on the details of 
the field being modeled. Some common examples are: 

(i) Point-mass field. If the Galaxy is represented as a 
point mass Mq at distance Ra, we have 

1 r.2 GMg . 

(ii) Isothermal field. For motion in a "halo" mass distri- 
bution modeled as an isothermal sphere (p ~ ^~^), with 
M(< r) = Mcir/Ra), we have 

a3 = -iai=n^ = ^. (B9) 
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(iii) Disk field. For motion in a disk described by local 
Oort constants A and B, with local density po, we have 



ai = -4A{A - B) 

a = A-B 



B^) 



(BIO) 

(Bll) 
(B12) 



In the (fairly good) approximation that the gravita- 
tional potential of the cluster stars may be represented close 
to the Jacobi surface simply as (t)c{r) ~ —GMtot/r, where 
Mtot is the cluster mass, the Jacobi radius may straightfor- 
wardly be shown to be 

--i^r"^ ,B,3, 

The ratio qs/qi determines the shape of the Jacobi surface. 
B1.5 Escaper Removal 

Stars are removed ( "stripped" ) from the system when they 
exceed a specified distance from the cluster center of mass 
(or density center) . For systems without an imposed Galac- 
tic tidal field, this stripping radius is arbitrary. For systems 
with a tidal field, the stripping radius is usually tied to the 
Jacobi radius of the cluster. For the runs described in this 
paper, stars were stripped when their distance from the clus- 
ter center exceeded twice the instantaneous Jacobi radius. 

B2 SeBa 

The stellar and binary evolution package SeBal^ is fully inte- 
grated into the kira integrator, although it can also be used 
as a stand-alone module for non-dynamical applications. 

B2. 1 Evolution of a single star 

Stars are evolved via the time dependent mass-radius rela- 
tions for solar metallicities given by Eggleton et al. (1989, 
with corrections by Eggleton et al. 1990 and Tout et 
al. 1997) . These equations give the radius of a star as 
a function of time and the star's initial mass (on the zero- 
age main-sequence — ZAMS). Neither the mass of the stellar 
core nor the rate of mass loss via a stellar wind are specified 
in this prescription. However, both quantities are impor- 
tant, both to binary evolution and to cluster dynamics. We 
include them using the prescriptions of Portegies Zwart & 
Verbunt (1996). 

Stars are subdivided within SeBa into the following 
types: 



brown dwarf 

main sequence 
Wolf-Rayet 

helium star 



planet Various types, such as gas gi- 

ants, etc.; also includes moons. 
Star 

with mass below the hydrogen- 
burning limit. 

Core hydrogen burning star. 

Massive (m > 25 M©) star which 

has lost its hydrogen envelope 

via a stellax wind. 
Helium core of a stripped giant, 

the result of mass transfer in a 
binary. Subdivided into helium 
core, carbon core and helium gi- 
ant. 

subgiant Hydrogen shell burning star. 

horizontal branch Helium core burning star. 

supergiant Double shell burning star. 

Thorne-Zytkow Shell burning hydrogen envelope 

with neutron star core, 
black hole Star with radius smaller than 

the event horizon. The result 

of evolution of massive (m > 
25 M0) star or collapsed neutron 
star. 

neutron star Subdivided into radio pulsar, X- 

ray pulsax and inert neutron star 
(m < 2M0). 

white dwarf Subdivided into helium dwarf, 

carbon dwarf and oxygen dwarf. 

disintegrated Result of Carbon detonation to 

Type la supernova. 

Stellar-wind mass loss is neglected for main-sequence 
stars with m < 25 M©. Following Langer (1998), more mas- 
sive stars lose mass with r'n oc m'^''' before becoming a Wolf- 
Rayet star (see Portegies Zwart et al. 1999, for the imple- 
mentation). These stars eventually collapse into black holes 
with mass rribh = 0.35mo — 12 Mq, where mo is the initial 
mass of the star. (For a star whose mass increases due to 
collisions or other processes, tuq is the highest mass reached 
by the star. The black hole radius equals the Schwarzschild 
radius: r = 2Gm/(? .) 

A star with a helium core mass between 2.2 and 5 M© 
becomes a neutron star. (These limits correspond to 8 M© 
and 25 M© ZAMS mass stars which evolve as isolated single 
stars.) At birth, a neutron star receives a velocity "kick" in 
a random direction. The magnitude of the velocity kick is 
chosen randomly from the distribution proposed by Hart- 
man (1997) 



P{u)du = — 



TT (l-hM2)2' 



(B14) 



II The name SeBa is taken from the ancient Egyptian word for 
'to teach', 'the door to knowledge' or '(multiple) star'. The exact 
meaning depends on the hieroglyphic spelling. 
** New equations which include mctallicity dependence have re- 
cently been made available by Hurley et al. (2000), and will be 
implemented in the next version of SeBa. 



with u = v/a and a = 600 kms ^. 

A star with a core mass less than 2.2 Mq sheds its 
envelope at the end of its evolution and becomes a white 
dwarf. The mass of the white dwarf equals the core mass of 
its progenitor at the tip of the asymptotic giant branch. 
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B2.2 Schematic evolution of a binary 

The evolution of a single isolated or unperturbed binary is 
carried out in the following steps (see sects. B2.3 and B2.4 
for details): 

The evolution of any binary stars by Determining the 
binary evolution timostop. This is the smallest timestep al- 
lowed by either of the stars. A stellar evolution timestep is 
1% of the time talcen for the star to evolve from the start of 
one evolutionary stage to the next — for example, from the 
zero-age main sequence to the terminal-age main sequence. 
(The stellar evolution step is not to exceed IGyear.) A list 
of these mile-posts along a star's evolutionary is provided 
in Sect. B2.1. A binary is evolved whenever one of its stars 
requires an update. 

If a binary is in a state of mass transfer (but not in a 
common envelope) the timestep is reduced such that < 1% 
of the donor' envelope is lost per step. 

1. Apply angular momentum loss by magnetic stellar 
wind. 

2. Apply angular momentum loss by gravitational wave 
radiation. 

3. Check for coalescence. 

4. Evolve primary star. 

(a) adjust binary parameters for stellar wind mass loss 

(b) resolve supernova. 

5. Check if binary still exists. The evolution of the pri- 
mary (or secondary) star may have resulted in a supernova 
which may disrupt the binary or resulted in a collision be- 
tween the two stars. 

6. Evolve secondary star. 

(a) Adjust binary parameters for stellar wind mass loss. 

(b) Resolve supernova. 

7. Check if binary still exists (see [5.]). 

8. Check for tidal circualrization and synchronization. 

9. Check if any star is Roche-lobe filling and identify the 
donor and the accretor. If no star is filling its Roche-lobe 
leave binary evolution and notify dynamics, otherwise pro- 
ceed with the following steps 

(a) Find moment mass transfer starts. 

(b) Check for binary stability 

if binary unstable apply commone envelope 

if components merge leave binary evolution and 
notify dynamics. 

(c) Calculate Cad, Cm and Cth. 

(d) Determine amount of mass loss from donor. 

(e) Determine amount of mass gained by accretor. 

(f) Subtrac mass from donor. 

(g) Add mass to accretor. Calculate new evolutionary 

state of accretor and Rejuvenate. 

(h) Calculate new binary parameters. 



B2.3 Evolution of binary parameters without mass 
transfer 

The orbital parameters of a binary are affected by the evo- 
lution of its components. We will not present here all the 
many details of binary evolution, but for clarity we sum- 
marize those which affect the dynamics or are important 
for interpreting our results. The details of the binary evolu- 
tion program SeBa are discussed in more detail by Portegies 
Zwart & Verbunt (1996) and Portegies Zwart & Yungelson 
(1998). 

Mass lost in a stellar wind is assumed to escape isotrop- 
ically from the binary system. If the companion accretes a 
fraction ^ of the other star's wind this implies 



a _ ^ Mo + nio 
fflo M + m 



(B15) 



Here Mo and mo are the initial primary and secondary mass, 
respectively, M and m are their final masses. And 



/ = 





r M 1 




( 


[mo\ 


mo / 



(B16) 



The fraction ^ is calculated via Bondi-Hoyle (1944) accre- 
tion assuming that the thermal velocity in the wind equals 
the escape velocity of the mass-losing star (for details see 
Portegies Zwart & Verbunt 1996). 

When the radius of one of the stars exceeds 5 times 
the orbital periastron separation a(l — e), orbital energy is 
transformed into oscillatory modes in the two stars. This 
leads to a decrease in the orbital separation and, due to 
conservation of angular momentum [a(l — e^) = constant], 
to the eventual circularization of the binary. 

Mass loss in a supernova is lost impulsively from the bi- 
nary system. As a result both the orbital separation and the 
eccentricity change: both increase if the pre-supernova or- 
bit was circular. The velocity kick (Eq. B14) received by the 
neutron star at formation is added randomly to its orbital 
velocity. New orbital parameters are then calculated assum- 
ing that the positions of the two stars are unchanged, mass 
is lost isotropically from the exploding star, and the com- 
panion is unaffected by the explosion. If the pre-supernova 
orbit is eccentric things get somewhat more complicated (see 
Portegies Zwart & Verbunt 1996). 

Low-mass stars may have magnetically coupled winds, 
and relatively large changes in angular momentum may oc- 
cur even though a negligible amount of mass escapes. We 
follow the prescription described by Rappaport et al. (1983) 
for tidally synchronized binaries in which at least one com- 
ponent is a main-sequence star or (sub) giant with mass 
0.7 < m/Mo < 1.5. 

Compact stars in short-period binaries and highly ec- 
centric binaries lose orbital energy and angular momentum 
via gravitational radiation. For such binaries, we use the ex- 
pressions provided by Peters (1964) to compute the time 
dependence of the orbital semi-major axis and eccentricity. 
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B2.4 Mass transfer in binaries 



Cad 



-0.221 - 2.847a; + 32.03x^ - 75.69a;* + 57.81x^(B21) 



When one star in a binary approaches its Roche- hniit, we 
iteratively detorniinc the moment at which contact occurs. 
The size of the Roche lobe is calculated as (Eggleton 1983) 

0.49 



r-Ri ■ 



0.6-|-g2/3in(l +g-i/3)' 



(B17) 



where q = m/M . The Rochc-lobo-fiUing star is then identi- 
fied as the donor and its companion as the accretor. 

B2.4.1. Unstable Mass transfer 

When a star fills its Roche lobe we first check for the possi- 
bility of Darwin- Riemann instability. This happens if 



Jdo 



(B18) 



where Jdonor and Jbin are the angular momenta of the Roche- 
lobe-filling star and the binary, respectively. 

During spiral-in the envelope of the donor is expelled, 
at the cost of orbital energy, following the prescription of 
Webbink (1984): 

^ Mc / l_f2ao Me 
ao ' ' 



a _ Mc C I + 2ao Mey 
ao M \ aX m J 



(B19) 



The parameters governing binary evolution are listed in Ta- 
ble B2. 

If the Roche-lobe-filling star is a main-sequence star 
or compact object, the two stars simply merge because the 
donor has no core-halo structure. If both the donor and the 
accretor are (sub)giants, we expel both envelopes at the cost 
of the binaries binding energy in a double inspiral. A merger 
occurs when the binary that remains after the common en- 
velope phase is semi-detached, in which case no more mass 
is lost (see sect. B4). 

B2.4.2. Stable mass transfer 

We calculate the time scale for mass transfer in a dynami- 
cally stable binary by considering the responses of the donor 
and the binary parameters to changes in the donor mass. For 
this purpose we define the logarithmic derivative 

^-(^),. <-«) 

for each of the following processes: 

Cad the change in donor radius due to adi-, 

abatic adjustment of hydrostatic equi- 
librium 

Cri the change in the size of the donor's, 

Roche lobe 

Cth the change in donor radius as it ad-, 

justs to a new thermal equilibrium 
The adopted values for Cad are as follows: For main se- 
quence stars with m > 0.7 we use Cad = 4 and half this 
value for lower-mass main sequence stars. For stars on the 
Hertzsprung gap and horizontal branch we use Cad = 2.25 
and 15, respectively. For other stars with a core-halo struc- 
ture (subgiants, supcrgiants and TZ objects), we use the 
following fit to the composite polytropic models of Hjellm- 
ing & Webbink (1987): 



where x = mcoTe/m. 

We use Cth between and 0.9 for main sequence stars 
and Cth = for all other stars, except those on the 
Hertzsprung gap and on the horizontal branch for which 
we use Cth = —2 and 15, respectively. 

The response of the Roche lobe to mass transfer Cri 
is calculated by transfering a infenitesimal amount of mass 
from the donor to the accreting star and study the response 
of the binary parameters. This test particle is transfered on 
the same timescale as was used in the previous mass transfer 
step. At first Roche-lobe contact, when there was no previ- 
ous mass transfer step, we assume that the test particle is 
transfered on a thermal timescale, which is a rather conser- 
vative choise. 

The time scale on which mass transfer proceeds is de- 
termined as follows: 

Cad < Cri dynamically unstable 

mass transfer proceeds on time 
scale Tdyn 

thermally unstable mass trans- 
fer proceeds on time scale Tth 
nuclear unstable mass trans- 
fer proceeds on time scale 

min(Tnuc,Tj), 

where the time scales associated with the various crite- 
ria are as follows: 



Cad > Cm and Cth < Cm 

Cad > Cri and Cth > Cri 



dynamic: 
thermal: 

nuclear: 

angular momentum loss: 



''"dyn 
7th 

rnuc 
TJ 



5.110-"A/r3/m [Myr] 
32m'^/{rL) [Myr] 

O.Ums 

*A>in/(^gr -t- </mb)- 



Here tms is a star's main-sequence lifetime, and m, r and L 
are its mass, radius and luminosity, respectively. The loss of 
angular momentum via gravitational radiation and magnetic 
braking are denoted by Jgr and Jmb, respectively. 

Table Bl gives a flavor of the various time scales on 
which mass transfer generally proceeds. However, the details 
depend critically on the orbital separation and on the mass 
and evolutionary state of both the donor and the accreting 
star. 

Some mass may be lost from the binary system during 
mass transfer. The new orbital parameters are calculated 
assuming that the mass lost from the binary carries specific 
angular momentum rjj (see Table B2). We calculate the final 
orbital separation using 



_o_ _ / Mm Y ( M + 
ao~ \ Momo ) V Mo -|- 



m \ 
mo) 



2'7j + l 



(B22) 



Here M = Mo — dM and m = mo + dm (so dM and dm are 
defined as positive quantities). The binary thus loses mass 

if dM — dm > 0. For the amount of mass accepted by the 
accretor, see Portegies Zwart & Verbunt (1996). 
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Table Bl. Schematic diagram of the time scales on which stable mass transfer (donor to accretor) proceeds. 
Donor: main sequence subgiant supergiant compact object 

Accretor: 

main sequence nucleajr/thermal thermal dynamic — 

(sub)giant — nuclear /thermal thermal/dynamic - 

compact object thermal/ami dynamic dynamic ami 



B3 Rejuvenation of the accretor 

An accreting star generally becomes more massive, which 

shortens its evolutionary time scale. The method described 
here is rather ad-hoc. We assume that a star accreting 5m 
of mass remains in the same evoutionary state (see the list 
in Sect. B2.1). The age of the star with mass m is t{m) and 
we want to know what is the age i(m + 5m) of the star 
with mass rn + 5m. At he moment the mass of the accretor 
increases from m to m + 5m the star is in evolutionary state 
i. It took the star ti{m) to reach that evolutionary state 
and this state lasts for Ti{m) = ti+i{m) — ti{m) for a star 
with mass m. The same stage for a star with mass m + 5m 
lasts for Ti{m-\-5m). The age of the star after accretion then 
becomes 

t{m+5m) = ti {m+5m) + 1^ ) ri(m+5m)7^.(B23) 

Here TZ is a fraction introduced to mimick the rejuvenation 
of the accretor {TZ > 1). The mass dumped on its surface 
may lead to some internal mixing, refreshing some of the 
helium core material with some of the freshly accreted Hy- 
drogen. This rejuvenation fraction is calculated with 



Table B2. Free parameters in binary evolution 



7^ 



m + (5m\ " 
m J 



(B24) 



and wc use TZ — 1 if the accreted material is not Hydrogen, 
in which case the accretor is not rejuvenated. Allowing TZ < 
1 would mimick that a star becomes older upon accreting 
material. The adopted value for k is listed in Table B2. 

We will give two examples of a m = 2 M© star which ac- 
cretes 5m = 0.2 Mq from a Helium rich companion {TZ = 1). 
For simplicity wc assume here that this amount of mass is 
transferred in an infcnitcsimal tiniestep. The main-sequence 
lifetime for a 2 M© star is about 801 Myeax and about 
608 Myear for a 2.2 M© star. If mass transfer starts at 
t = 700 Myear, the 2 M© accretor is still on the main se- 
quence. After mass transfer the accreting star has an age of 
531 Myear and is still on the main sequence. 

If mass transfer started at t = 1 Gyear things become 
somewhat more complicated. The 2 M© accretor is then on 
the horizontal branch. The time it takes from zero-age to 
the beginning of the horizontal branch is about 938 Myear, 
for a 2.2 M© star this is about 712 Myear. The 2M0 star 
spends roughtly 84 Myear on the horizontal branch, where a 
2.2 M© star spends only 82 Myear in that stage. Substitution 
of these numbers into Eq. B23 results in an age of the post 
mciss transfer star of 773 Myear. 



term 


value 


description 


K 


1 


accretor rejuvenation factor 


VJ 


2 


specific angular momentum loss per unit mass 


X 


0.5 


envelope binding energy fraction 


ace 


4 


common envelope constant 



B4 Result of a merger or collision 

We have adopted a set of simple prescriptions to specify 
the outcome of stellar collisions. In the future these pre- 
scriptions can be refined when more accurate calculations 
become available. As a rule of the thumb, the result of a 
collision is the conservative accretion of the lower-mass star 
onto the more massive star. The accretor will then be reju- 
venated as described in Eq. B24. This rule is violated when 
one component is a giant or a compact object. A detailed 
prescription of how to calculate the evolutionary state of 
such a merger is presented in paper I. 

We describe our treatment of the possible outcomes of 
encounters between two stars, ordered by the evolutionary 
state of the more massive of the two (the primary). Table B3 
summarizes this treatment. 



34.1 Main-sequence primary 

If both stars involved in the encounter are main-sequence 

stars, then the less massive star is accreted conservatively 
onto the more massive star. The resulting star is a reju- 
venated main-sequence star (see Lai et al. 1993, Lombardi 
et al. 1995). The details of this procedure are described in 
Appendix C4 of Portegies Zwart & Verbunt (1996). 

If the less massive star in the encounter has a well de- 
veloped core (giant or subgiant), this core becomes the core 
of the merger product. The main-sequence star and the en- 
velope of the giant are combined to form the new envelope 
of the merger. In general, the mass of the core is relatively 
small compared to the mass of the envelope, and the star is 
assumed to continue its evolution through the Hertzsprung 
gap. Note that this type of encounter can only occur when 
the main-sequence star is itself a collision product (e.g. a 
blue straggler). 

When a main-sequence star encounters a less massive 
white dwarf, we assume that the merger product is a gi- 
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Table B3. Simplified representation of possible merger outcomes. 
The four columns correspond to the four choices given for the 
type of massive star (primary), while the four rows indicate the 
type of less massive star (secondary): main-sequence star (ms), 
(sub)giant (sg), white dwarf (wd) and neutron star (ns). In this 
table we do not distinguish between stars in the Hertzsprung gap 
(Hg) or on the first and second ascent on the asymptotic-giant 
branch (AGB). 

primary 



star 


ms 


sg 


wd 


ns 


ms 


ms 


sg 


wd 
+ 
disc 


IIS 

+ 
disc 


sg 


Hg 


AGB 


wd 
+ 
disc 


ns 
+ 
disc 


wd 


sg 


AGB 







TZO TZO 



age 11.34 Gyr with a 0.45 M© main- sequence star produces 
a giant of 1.4 Mq with an age of 2.67 Gyr. 

When both stars are (sub)giants, the two cores axe 
merged and form the core of the merger product (see Davies 
ct al. 1991 and Rasio & Shapiro 1995). Half the envelope 
mass of the less massive star is accreted onto the primary. 
The merger product continues its evolution starting at the 
next evolutionary state — a (sub)giant continues its evolu- 
tion on the horizontal branch, and a horizontal branch star 
becomes an asymptotic-giant branch star. The reasoning be- 
hind this assumption is that an increased core mass corre- 
sponds to a later evolutionary stage. 

If the less massive star is a white dwarf, then its mass 
is simply added to the core mass of the giant, and the enve- 
lope is retained. If the age of the giant before the encounter 
exceeds the total lifetime of a single unperturbed star with 
the mass of the merger, then the newly formed giant im- 
mediately sheds its onvolopo and its core turns into a single 
white dwarf. Otherwise the merged giant is assumed to have 
the same age (in years) as the giant before the collision, and 
continues its evolution as a single unperturbed star. 

If the other star is a less massive neutron star, a Thorne- 
Zytkow object is formed. 



ant whose core and envelope have the masses of the white 
dwarf and the main-sequence star, respectively. We then de- 
termine its evolutionary state as follows. We calculate the 
total time tagb that a single, unperturbed star with mass 
equal to that of the merged star would spend on the asymp- 
totic giant-branch, and the mass mc,agb of its core at the 
tip of the giant branch. The age of the merger product is 
then calculated by adding tagb"^c/?«c,agb to the age of an 
unperturbed star with the same mass at the bottom of the 
asymptotic giant branch. For example, a single, unperturbed 
1.4Mq star leaves the main-sequence after 2.52 Gyr, spends 
60 Myr in the Hertzsprung gap, moves to the horizontal 
branch at 2.96 Gyr, and reaches the tip of the asymptotic 
giant branch after 3.06 Gyr, with a core of 0.64 Mq. Thus, if 
a 0.6 M0 white dwarf merges with an 0.8 Mq main-sequence 
star, the merger product has an age of 2.87 Gyr, leaving it 
another 180 Myr before it reaches the tip of the asymptotic 
giant-branch. 

If the less massive star is a neutron star or black hole a 
Thorne-Zytkow object (1977) is formed. 



B4.2 Evolved primary 

When a (sub) giant or asymptotic branch giant encounters 
a less massive main-sequence star, the main-sequence star 
is combined with the envelope of the giant, which stays 
in the same evolutionary state. Its age within that state 
is changed, however, according to the rejuvenation calcu- 
lation described in Sect. C3 of Portegies Zwart & Verbunt 
(1996). For example, an encounter of a giant of 0.95 M0 and 



B4-3 White-dwarf primary 

In an encounter between a white dwarf and a less massive 
main-sequence star, the latter is assumed to be completely 
disrupted, and forms a disk around the white dwarf (Rasio 
& Shapiro 1991). The white dwarf accretes from this disk 
at a rate of 1 percent of the Eddington limit. If the mass 
in the disc exceeds 5% of the mass of the white dwarf, the 
excess mass is expelled from the disc at a rate equal to the 
Eddington limit. 

If a white dwarf encounters a less massive (sub)giant, 
a new white dwarf is formed with a mass equal to the sum 
of the prc-encountcr core of the (sub)giant and the white 
dwarf. The newly formed white dwarf is surrounded by a 
disk formed from half the envelope of the (sub)giant before 
the encounter. The factor half is rather arbitrary and based 
on the lack of detailed calculations which provide a proper 
rmmber. If the mass of the white dwarf exceeds the Chan- 
drasekhar limit it explodes in a type la supernova, leaving 
no remnant (Nomoto & Kondo 1991; Livio & Truran 1995). 

A collision between two white dwarfs results in a sin- 
gle white dwarf with mass equal to the sum of the original 
masses. If the total mass of the collision product exceeds the 
Chandrasekhar mass, it explodes in a type la supernova. 

Collisions between white dwarfs and neutron stars or 
black holes result in the formation of an accretion disc 
around the compact object; the white dwarf is destroyed. 
Following the accretion, the neutron star may collapse into 
a black hole. 
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B^.^ Neutron-star or black-hole primary 

All collisions involving a neutron star or black hole primary 
lead to the formation of a massive disk around the compact 
star. If the compact star had a disk prior to the collision, 
this disk is expelled. This disk accretes onto the compact 
star. Wo chose, rather arbitrarily that the accretion rate is 
5% of the Eddington limit. An accreting neutron star turns 
into a millisecond radio pulsar, or — when its mass exceeds 
2M0 — into a black hole. 

B5 Communication between SeBa and kira 

Due to the interaction between stellar evolution and stellar 
dynamics, it is difficult to solve for the evolution of both 
systems in a completely self-consistent way. The trajectories 
of stars are computed using a block timestep scheme, as 
described earlier. Stellar and binary evolution is updated at 
fixed intervals (every 1/64 of a crossing time, typically a few 
thousand years) . Any feedback between the two systems may 
thus experience a delay of at most one timestep. Internal 
evolution time steps may differ for each star and binaxy, and 
depend on binary period, perturbations due to neighbors, 
and the evolutionary state of the star. Time steps in this 
treatment vary from several milliseconds up to (at most) a 
million years. 

After each 1/64 of a crossing time, all stars and bina- 
ries are checked to determine if evolutionary updates are 
required. Single stars are updated every 1/100 of an evo- 
lution timestep or when the mass of the star has changed 
by more than 1% since the last update. A stellar evolution 
timestep is the time taken for the star to evolve from the 
start of one evolutionary stage to the next (sec sect. B2.2). 

After each stellar evolution step the dynamics is no- 
tified of changes in stellar radii, but changes in mass are, 
for reasons of efficiency, not passed back immediately (mass 
changes generally entail recomputing the accelerations of 
all stars in the system). Instead, the "dynamical" masses 
are modified only when the mass of any star has changed 
by more than 1%, or if the orbital parameters, semi-major 
axis, eccentricity, total mass or mass ratio of any binary has 
changed by more than 0.1%. 

B6 Mass loss from stars and binaries 

Fast (sudden) and slow (gradual) mass loss affects the dy- 
namics of the stellar system in different ways. Mass loss is 
considered fast when it takes place within a fraction of an 
orbital time scale. For single stars, this time scale is on the 
order of the crossing time of the star cluster. For binaries, 
it is much shorter — on the order of the binary orbital pe- 
riod. Mass loss during a supernova explosion is considered 
fast, stellar winds and mass lost from a binary during mass 
transfer are considered slow. 

Due to the discretized time steps of the stellar dynamics 
and the stellar evolution, from the point of view of the dy- 
namics mass is lost in "bursts." For example, an asymptotic 



giant star with a strong stellar wind may lose its entire enve- 
lope in a hundred steps spanning roughly one crossing time, 
while a supergiant might lose its entire envelope instanta- 
neously in a supernova. Mass loss for single stars affects the 
dynamics of the entire stellar system. For binaries and mul- 
tiple systems, mass loss from a member star directly affects 
the orbital characteristics of its neighbors. 

The rate of mass loss is particulaxly important for bi- 
naries. Slow mass loss via a stellar wind will soften a binary 
system, but will not affect its eccentricity or its center of 
mass velocity. (This is true if the binary is unperturbed. In 
a perturbed binary, the eccentricity and center of mass ve- 
locity are both affected by stellar wind mass loss.) Sudden 
mass loss, on the other hand, can dramatically affect the 
binary's internal parameters. For unperturbed binaries, the 
effects of mass loss from both component stars are computed 
consistently using SeBa. Changes in binary parameters are 
calculated and the dynamics is notified, thereby transmit- 
ting the information to the rest of the stellar system via the 
integrator. 

For perturbed binaries and multiples (and also hierar- 
chical systems where the inner binary is unperturbed), the 
integrator takes care of the dynamical effects of stellar mass 
loss. By construction, mass transfer cannot occur in a per- 
turbed binary or multiple system. If a supernova occurs in a 
perturbed binary, any slow mass loss is accounted for before 
fast mass loss occurs, since a star which is about to explode 
generally loses a significant fraction of its mass in a stellar 
wind before the supernova event itself Supernova remnants 
do not lose mass. This assumption breaks down when the 
binary companion of the exploding stax loses a significant 
fraction of its mass between the moment of the supernova 
and the end of the stellar update timestep. (This can hap- 
pen if the binaxy companion is either a Wolf-Rayet star or a 
supergiant.) The stellar evolution time steps in these cases 
are taken sufficiently small (on the order of a hundred years) 
to ensure that this causes a negligible error. 



B7 Collisions and mergers 

We draw a distinction between "mergers" and "collisions." A 
merger may result from mass transfer or a common-envelope 
phase during the evolution of an unperturbed binaxy. The 
binary node is then replaced by the merger product. The 
product of a merger is generally different from the result of 
a collision, since a merger is often preceded by a phase of 
mass transfer which affects the masses of both stars. 

Collisions may occur between between single stars 
(which are part of a binary tree) or between stars in a per- 
turbed binary. Since the integrator may miss the precise mo- 
ment of closest approach, the orbital elements of each "close" 
pair of stars is calculated after each integration step. A colli- 
sion occurs when the staxs axe found to be within dcoii times 
the sum of their radii at periastron: p < dcoii(ri -I- r2). In 
this case, the two stars are replaced by the collision product, 
which is placed in the center of mgiss and with the center of 
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mass velocity of the original two-body system. The charac- 
teristics of the collision product are calculated using SeBa 
(Sect.B4). 

A collision may also occur when an unperturbed binary 
in a state of mass transfer is perturbed by a close encounter 
with another cluster member. Such a induced collision may 
be triggered by a close flyby, or in a multiple system with a 
perturbed outer orbit. The collision occurs if the sum of the 
component radii exceeds the distance between the two stars 
at the moment the binary becomes perturbed. 
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